Samantha Brunel

E Andrew

Andrew Bennett

Laurent Cardin

Damien Garraud

Hélène Barrand Emam

Alexandre Beylier

Bruno Boulestin

Fanny Chenal

Elsa Ciesielski

Fabien Convertini

E Andrew Bennett

Bernard Dedet

Stéphanie Desbrosse-Degobertiere

Sophie Desenne

Jerôme Dubouloz

Henri Duday

Gilles Escalon

Véronique Fabre

Eric Gailledrat

Muriel Gandelin

Yves Gleize

Sébastien Goepfert

Jean Guilaine

Lamys Hachem

Michael Ilett

François Lambach

Florent Maziere

Bertrand Perrin

Suzanne Plouin

Estelle Pinard

Ivan Praud

Isabelle Richard

Vincent Riquier

Réjane Roure

Benoit Sendra

Corinne Thevenet

Sandrine Thiol

Elisabeth Vauquelin

Luc Vergnaud

Thierry Grange
email: thierry.grange@ijm.fr

Eva-Maria Geigl
email: eva-maria.geigl@ijm.fr

Melanie Pruvost
email: melanie.pruvost@u-bordeaux.fr

Ancient genomes from present-day France unveil 7,000 years of its demographic history

Keywords: paleogenomics, migration, Neolithic, Mesolithic, Protohistory, population genomics

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

Introduction

Over the last 10,000 years, populations in Western Eurasia have undergone two major cultural shifts: from a hunter-gatherer lifestyle transitioning to a lifestyle based on food production during the Neolithic [START_REF] Bellwood | First Farmers: the Origins of Agricultural Societies[END_REF] , and the development and improvement of metallurgy over the course of the third and second millennia BCE, giving rise to the Bronze Age that in turn evolved into the Iron Age during the last seven centuries BCE [START_REF] Kristiansen | The Rise of Bronze Age Society: Travels, Transmissions and Transformations[END_REF] .

Ancient genomes have been instrumental in characterizing past populations [START_REF] Fu | The genetic history of Ice Age Europe[END_REF][START_REF] Keller | New insights into the Tyrolean Iceman's origin and phenotype as inferred by whole-genome sequencing[END_REF][START_REF] Lazaridis | Ancient human genomes suggest three ancestral populations for present-day Europeans[END_REF][START_REF] Olalde | Derived immune and ancestral pigmentation alleles in a 7,000-year-old Mesolithic European[END_REF][START_REF] Sánchez-Quinto | Genomic affinities of two 7,000-year-old Iberian huntergatherers[END_REF][START_REF] Skoglund | Origins and Genetic Legacy of Neolithic Farmers and Hunter-Gatherers in Europe[END_REF] . Diachronic series have investigated temporal genome-wide dynamics within defined European regions such as Central Europe [START_REF] Gamba | Genome flux and stasis in a five millennium transect of European prehistory[END_REF] , Great Britain 10 or Iberia [START_REF] Villalba-Mouco | Survival of Late Pleistocene Hunter-Gatherer Ancestry in the Iberian Peninsula[END_REF][START_REF] Valdiosera | Four millennia of Iberian biomolecular prehistory illustrate the impact of prehistoric migrations at the far end of Eurasia[END_REF] , revealing that these cultural changes were the result of demographic changes that deeply altered the genetic makeup of Western Eurasian populations [START_REF] Lazaridis | Ancient human genomes suggest three ancestral populations for present-day Europeans[END_REF][START_REF] Gamba | Genome flux and stasis in a five millennium transect of European prehistory[END_REF][START_REF] Olalde | The Beaker phenomenon and the genomic transformation of northwest Europe[END_REF][START_REF] Haak | Massive migration from the steppe was a source for Indo-European languages in Europe[END_REF] .

In France, the demographic processes underlying these transitions are not yet explored at a territory-wide scale. Only a handful of studies restricted to individual archeological sites have been performed, relying on partial mitochondrial sequence information [START_REF] Beau | Multi-scale ancient DNA analyses confirm the western origin of Michelsberg farmers and document probable practices of human sacrifice[END_REF][START_REF] Lacan | Ancient DNA suggests the leading role played by men in the Neolithic dissemination[END_REF][START_REF] Rivollat | When the Waves of European Neolithization Met: First Paleogenetic Evidence from Early Farmers in the Southern Paris Basin[END_REF][START_REF] Rivollat | Ancient mitochondrial DNA from the middle neolithic necropolis of Obernai extends the genetic influence of the LBK to west of the Rhine[END_REF] or on partial Y chromosome sequences [START_REF] Lacan | Ancient DNA suggests the leading role played by men in the Neolithic dissemination[END_REF] . Due to its geographic location, France occupies a strategic position in the understanding of population migrations in Western Europe [START_REF]La protohistoire de la France[END_REF][START_REF] Demoule | La révolution néolithique en France. The neolithic revolution in France[END_REF][START_REF] Constantin | Anneaux d'origine méridionale dans le Rubané récent du Bassin parisien. Rings of southern origin in the Recent Rubane of the Paris Basin[END_REF] . Both the Central European and Mediterranean currents of neolithization (the Linearbandkeramik (LBK) complex in the North and the Impressa and Cardial complexes in the South) became established there, although the extent of their interactions remains an open question [START_REF]La protohistoire de la France[END_REF][START_REF] Demoule | La révolution néolithique en France. The neolithic revolution in France[END_REF][START_REF] Constantin | Anneaux d'origine méridionale dans le Rubané récent du Bassin parisien. Rings of southern origin in the Recent Rubane of the Paris Basin[END_REF] . Moreover, the consequences on the gene pool of the region during the development and spread of metallurgy across Europe, at the beginning of the Bronze and Iron Ages, are still unknown.

To investigate the demographic dynamics over the 7,000-year-transect from the Mesolithic period, before the onset of agriculture, to the Iron Age, we genetically analyzed 243 unique individuals sampled from 54 different archeological sites, enriching for both complete mitochondrial genomes and a panel of 120 nuclear single nucleotide polymorphisms (SNPs) including Y chromosome SNPs, as well as a partially overlapping dataset of 58 low coverage genomes. Our data from ancient Western Europe fill a geographic gap that was so far missing, allowing a more global view of past population dynamics in Europe.

Results and Discussion

All analyzed individuals belong to well-defined archaeological contexts from the Mesolithic period to the Iron Age and 40 were directly radiocarbon dated for this study (see SI Appendix, Table S1-2). We built sequencing libraries for 243 individuals and enriched them for mitochondrial genomes and a panel of 120 selected SNPs across the nuclear genome (see SI Appendix, Table S1-1). Along with complete mitochondrial genomes of 223 individuals and Y chromosome haplogroups of 62 male individuals, we report 58 newly generated ~0.05X to 0.5X genomes sequenced from UDG-treated libraries. We combined our genotype data with either complete genomes or 1240k genome-wide polymorphisms in ancient individuals, and genotypes of modern individuals from Europe, Caucasus and the Near East genotyped on the Affymetrix Human Origins array [START_REF] Fu | The genetic history of Ice Age Europe[END_REF][START_REF] Lazaridis | Ancient human genomes suggest three ancestral populations for present-day Europeans[END_REF][START_REF] Gamba | Genome flux and stasis in a five millennium transect of European prehistory[END_REF][START_REF] Olalde | The Beaker phenomenon and the genomic transformation of northwest Europe[END_REF][START_REF] Valdiosera | Four millennia of Iberian biomolecular prehistory illustrate the impact of prehistoric migrations at the far end of Eurasia[END_REF][START_REF] Haak | Massive migration from the steppe was a source for Indo-European languages in Europe[END_REF][START_REF] Lipson | Parallel palaeogenomic transects reveal complex genetic history of early European farmers[END_REF][START_REF] Mathieson | Genome-wide patterns of selection in 230 ancient Eurasians[END_REF] . Complete mitochondrial genomes were retrieved from this ancient dataset and were grouped based on temporal and geographic proximity, often translating into cultural classifications (see SI Appendix S3).

Contamination estimates derived from X chromosome polymorphisms in male individuals were consistently low (between 0.68%±1.26% and 3.55%±0.62% depending on the method considered), with only one individual displaying more than 5% X chromosome contamination (PIR3037AB) (see SI Appendix Table S1-1).

The French Mesolithic substrate

We generated genomic data for five Mesolithic individuals from Les Perrats cave (Agris, Charente) dated from 7177-7057 cal. BCE including low coverage shotgun sequencing for three of them (0.111-0.134X). All possess haplogroup U5b, demonstrating in the territory of present-day France the presence of a Mesolithic substrate akin to that described elsewhere where it was characterized by haplogroups U5b (85%) and, to a lesser extent, U5a (15%) [START_REF] Fu | The genetic history of Ice Age Europe[END_REF][START_REF] Posth | Pleistocene Mitochondrial Genomes Suggest a Single Major Dispersal of Non-Africans and a Late Glacial Population Turnover in Europe[END_REF] .

Mesolithic hunter-gatherers discovered in France are located near the edge of a PCA plot obtained from genome-wide genotypes, falling near the Mesolithic individuals from Western Europe [START_REF] Lazaridis | Ancient human genomes suggest three ancestral populations for present-day Europeans[END_REF][START_REF] Olalde | Derived immune and ancestral pigmentation alleles in a 7,000-year-old Mesolithic European[END_REF][START_REF] Gamba | Genome flux and stasis in a five millennium transect of European prehistory[END_REF] (Figure 1B and Figure S4-3). While the ancestry from a ~14,000-year-old individual from Villabruna in Italy was found to be dominant in Holocene hunter-gatherers across Western and Central Europe

Successive waves of migration and admixture over the course of the Neolithic in France

The arrival of an Anatolian Neolithic-associated ancestry component at the onset of the Neolithic is clearly visible in our dataset from both uniparentally-inherited markers and at the genome-wide level.

Maternal lineages from Neolithic farmers in France (from ~5300 BCE onwards) are more diverse than those from HGs, and display variable frequencies across the time transect (see SI Appendix, Figure S3-2). Early and Middle Neolithic individuals share more affinity with present-day Southern Europeans (see SI Appendix, Figure S4-2), and their genetic variation is encompassed within that of contemporaneous European populations (Figure 1B). Three of our individuals associated with the LBK culture, which stems from the Danubian Neolithic current, cluster genetically with other Central European Early Neolithic (EN) individuals and share more drift with other early farmers associated with the LBK cultural complex (see SI Appendix, Figure S4-6). Values of D (Mbuti, France_EN;

Western Hunter-Gatherers (WHG), Anatolia_Neolithic) are consistent with those observed in LBKassociated individuals, with the exception of Mor6 (~7,100 calBP) in north-eastern France (Figure 2B). Despite a consistent dating and cultural assignment to the LBK, this individual falls within the genetic diversity of Iberian_EN individuals and harbors the highest proportion of shared alleles with WHG. After testing the potential origin of this ancestry with qpAdm, we found the best fit to be with a two-source ancestry model between Anatolia_Neolithic and GoyetQ2 (p-value: 0,0911155), a situation so far only described in Early Neolithic individuals from the Iberian Peninsula. This is the northern-most description of such ancestry. To date, the distribution of the Magdalenian-associated ancestry in the hunter-gatherer population in Western Europe that admixed with the Neolithic farmers in Europe is not clear because of the scarcity of genetic data from Mesolithic HGs. The observation of this ancestry in the Iberian Peninsula has been interpreted as a feature distinguishing the Cardial from the LBK Neolithic migrants [START_REF] Valdiosera | Four millennia of Iberian biomolecular prehistory illustrate the impact of prehistoric migrations at the far end of Eurasia[END_REF] .

The observation that the ancestry of GoyetQ2 was found alongside Villabruna ancestry in HGs from Les Perrats in western France (~9,100 calBP), as well as in HGs from Rigney1 (central-eastern France;

~15,500 cal. BP), and from two caves in south-western Germany, Hohlefels (~15,000 calBP) and Burkhardtshöhle (~14,600 calBP), rather indicates that this mixed ancestry is characteristic of Western European HGs [START_REF] Villalba-Mouco | Survival of Late Pleistocene Hunter-Gatherer Ancestry in the Iberian Peninsula[END_REF] . Thus, the admixture between HGs with GoyetQ2 ancestry and Neolithic migrants is not necessarily characteristic of the Cardial Neolithic migration wave as it seems to have occurred as well in more northern locations in Western Europe.

While the only male Mesolithic individual from our dataset was assigned to Y haplogroup I, LBK individuals from France showed a relative diversity, belonging to Y haplogroups C1a2, G2a and H2. GoyetQ2-like ancestry, suggesting admixture with different populations of hunter-gatherers (Figure 1C).

The Bell Beaker period in France

At the end of the Neolithic, the Bell Beaker complex (named after the particular shape of its ceramics) appears and becomes widespread across Europe. The Beaker phenomenon is peculiar because of its large geographic distribution and because of its coexistence with local Late Neolithic and Copper Age cultures.

We report two Bell Beaker-associated individuals (CBV95 and PEI2), that we co-analyzed with previously reported contemporaneous individuals from Europe, including France Anatolia_Neolithic, Villabruna and Yamnaya_Samara, we could, however, detect 28.3% of Steppe ancestry in PEI2 (Figure 1C). These observations are consistent with previous findings and confirm that Steppe ancestry appeared later and with a lower impact in Southwestern Europe than in other parts of the continent [START_REF] Olalde | The Beaker phenomenon and the genomic transformation of northwest Europe[END_REF][START_REF] Olalde | The genomic history of the Iberian Peninsula over the past 8000 years[END_REF] . Surprisingly, the admixture model did not include Villabruna as an ancestry source for either CBV95 and PEI2, which differs from previously known Late Neolithic individuals.

Relative continuity between the Bronze Age and the Iron Age

Bronze Age individuals in France display new mitochondrial haplogroups such as U2, U4 and [START_REF] Olalde | The Beaker phenomenon and the genomic transformation of northwest Europe[END_REF][START_REF] Olalde | The genomic history of the Iberian Peninsula over the past 8000 years[END_REF][START_REF] Allentoft | Population genomics of Bronze Age Eurasia[END_REF] .

Analysis of genetic markers associated to phenotypes

To follow the evolutionary timelines of genetic adaptation to changing lifestyles and environments, we analyzed 73 autosomal loci, associated with both physical and physiological traits, genotyped across 149 individuals (Figure 3, see SI Appendix Table S2-1 for details on the targeted positions). Among the phenotypic traits of interest, we focused on genetic variants involved in eye and skin pigmentation (such as SLC24A5, SLC45A2, GRM5, HERC2, IRF4, TYR). It has been previously

shown that Mesolithic HGs and Anatolian Neolithic farmers differ in their pigmentation alleles, with the latter carrying derived alleles responsible for a light skin that are near fixation in Europe today [START_REF] Mathieson | Genome-wide patterns of selection in 230 ancient Eurasians[END_REF] , while HGs were described as having a rather dark skin and light eye color [START_REF] Olalde | Derived immune and ancestral pigmentation alleles in a 7,000-year-old Mesolithic European[END_REF][START_REF] Brace | Population Replacement in Early Neolithic Britain[END_REF] . The Mesolithic individuals from whom we retrieved reliable genotype information after stringent filtering (1 to 3 individuals) possessed the ancestral pigmentation variant for loci SLC45A2 and GRM5, which had been associated with darker skin. In contrast, the derived alleles of these two loci are found at significantly higher frequency in the French Neolithic population (SLC45A2: 24.9%, GRM5: 38.4%), although they remain far from the frequency reached by these alleles in present-day Europe, where they are near fixation in the population (93.8 and 68.9%, respectively). SLC24A5, described as the principal mutation associated with light skin pigmentation, is found in 96% in our Neolithic population, close to the 99% found in present-day Europe. These results suggest different evolutionary timelines for these mutations in France, where a similar phenotype (light skin pigmentation) possibly resulted from several selection events, although the genetic determinism of these traits is likely more complex and not encoded by these three loci alone [START_REF] Hysi | Genome-wide association meta-analysis of individuals of European ancestry identifies new loci explaining a substantial fraction of hair color variation and heritability[END_REF] . The Mesolithic individual that displayed coverage of variants involved in the pigmentation of the eye color (HERC2 and IRF4), carried the alleles associated with the blue eye color. At the time of the Neolithic, the frequency of these two variants reaches 29.6% and 24.4%, respectively. Interestingly, these two variants seemingly have different trajectories after the Neolithic. While HERC2 rises in frequency again from 37.5% during the Bronze Age to finally reach more than 63% in present-day Europe, IRF4 further decreases from 23.81% during the Bronze Age to reach a frequency of 11% at present. None of the genotyped Neolithic individuals carried the mutation responsible for the persistence of lactase in Europeans, consistent with a later origin for this mutation [START_REF] Haak | Massive migration from the steppe was a source for Indo-European languages in Europe[END_REF][START_REF] Allentoft | Population genomics of Bronze Age Eurasia[END_REF] . We also analyzed several innate immunity-related variants that show signs of recent positive selection in modern European populations [START_REF] Olalde | Derived immune and ancestral pigmentation alleles in a 7,000-year-old Mesolithic European[END_REF][START_REF] Hunt | Newly identified genetic risk variants for celiac disease related to the immune response[END_REF][START_REF] Yasuda | Association of single nucleotide polymorphisms in endothelin family genes with the progression of atherosclerosis in patients with essential hypertension[END_REF][START_REF] Zhou | Meta-analysis of genetic association of chromosome 9p21 with early-onset coronary artery disease[END_REF] . For a number of genes involved in the immune response (NOS2A, TOLLIP, CCL18, STAT6, IL3), the frequency of derived alleles was comparable to that found in present-day Europeans, indicating that selection on these loci predates the Neolithic, whereas frequencies of variants in the TLR1-6-10 locus, the Toll-Like receptor family involved in the innate immune response, are much lower than those found in present-day Europeans [START_REF] Barreiro | Evolutionary dynamics of human Toll-like receptors and their different contributions to host defense[END_REF] . The TLR1-6-10 gene cluster is possibly associated with resistance to leprosy, tuberculosis or other mycobacteria [START_REF] Barreiro | Evolutionary dynamics of human Toll-like receptors and their different contributions to host defense[END_REF] . Our data suggest that the onset of agriculture and the proximity with domesticated animals did not exert a selective pressure on this locus. Likewise, variants associated with coeliac disease in SH2B3 or SLC22A4 that could have been the target for positive selection owing to the role they play in inflammations and the elimination of a wide array of environmental toxins [START_REF] Hunt | Newly identified genetic risk variants for celiac disease related to the immune response[END_REF] , did not reach yet the frequency of modern populations.

Conclusion

The study of 7000 years of demographic history through the study of ancient genomes throughout France reveals successive migration events mostly accompanied by major cultural changes and admixture events whose traces are still found in current European populations.

Our genomic data obtained for three Mesolithic individuals from Western France show that the late survival of the Upper Paleolithic Magdalenian-associated ancestry (represented by GoyetQ2)

was not restricted to the Iberian Peninsula. This finding expands the region where this ancestry is found and raises the question of where this admixture occurred. More Mesolithic individuals from Western Europe need to be studied to better characterize the history of this ancestry.

While the first Neolithic people who inhabited the territory of present-day northern France descended from early farmers in Anatolia, therefore forming a group distinct from the local huntergatherers, we could detect later admixture with hunter-gatherers supporting their progressive integration within farming communities. Although widespread, the process of admixture with huntergatherers shows regional variability and origins. However, the scarcity of the data available for the South of France does not allow us to conclude whether contact occurred between the LBK and Cardial complexes, nor to identify any legacy that may be present in the Neolithic cultures that followed.

Genomic data from other French regions in the North and the South will be necessary to confirm the hypotheses on the origin and the relations of the different cultural groups of the Middle and Late

Neolithic.

The study of a selection of nuclear markers revealed differences between Neolithic people and present-day Europeans. Derived allele frequencies indicate selection on loci involved in pigmentation, diet and immunity that are consistent with adaptation to high latitudes and changes in diet. Finally, neither Mesolithic HGs nor Neolithic farmers on the territory of present-day France carried the allele allowing for lactase persistence during adulthood, in agreement with studies of Neolithic populations further East that reported the presence of this allele only by the end of the Bronze Age [START_REF] Haak | Massive migration from the steppe was a source for Indo-European languages in Europe[END_REF][START_REF] Allentoft | Population genomics of Bronze Age Eurasia[END_REF] .

An outcome of our study is the genetic continuity between Bronze and Iron Age individuals of France but with less heterogeneity between individuals from Iron Age. In such context, tracing migrations from one part of Europe to the other becomes more challenging and will require a deeper temporal and geographical resolution.

Material and methods

Processing of archaeological samples

A total of 243 individuals from 54 archeological sites across France passed screening for the preservation of ancient DNA and were processed for sequencing (see SI Appendix Table S1, Figure S1-1). These individuals covered a time period ranging from the Mesolithic to the Iron Age, with the majority of samples covering the Neolithic period (~5500-3500 BCE). As direct dates were not available for all samples, they were labeled and grouped for analyses based on geographical proximity and cultural assignment (see SI Appendix Figure S1-1).

Targeted DNA enrichment and sequencing DNA was extracted from petrous bones and teeth and used to build double-stranded Illumina libraries with thorough decontamination procedures and optimized protocols in a dedicated facility at the Institut Jacques Monod (see SI Appendix SI.2) [START_REF] Massilani | Past climate changes, population dynamics and the origin of Bison in Europe[END_REF] . Sequence enrichment for both the whole mitochondrial genome and nuclear regions was performed through a capture approach using biotinylated RNA baits obtained through in vitro transcription of PCR amplified target regions ( 37 , see SI Appendix SI.2 for details on the procedure). Both enriched and whole genome shotgun libraries were sequenced on an Illumina MiSeq or Illumina NextSeq at the Institut Jacques Monod in Paris, and the Institut de Recherche Biomédicale des Armées in Brétigny-sur-Orge, respectively.

Data processing and authentication

Adapters were trimmed and overlapping paired-end reads were merged using LeeHom [START_REF] Renaud | leeHom: adaptor trimming and merging for Illumina sequencing reads[END_REF] .

Merged reads were then quality trimmed with cutadapt and only reads longer than 30bp were retained.

Mapping was performed with BWA version 1.2.3 [START_REF] Villalba-Mouco | Survival of Late Pleistocene Hunter-Gatherer Ancestry in the Iberian Peninsula[END_REF] against the Cambridge reference sequence (rCRS) [START_REF] Andrews | Reanalysis and revision of the Cambridge reference sequence for human mitochondrial DNA[END_REF] for human mitochondrial DNA with a duplication of its first 100 bases at the end to ensure mapping of the reads overlapping the junction resulting from the virtual linearization of the circular mitogenome, and to the human reference genome hs37d5. In both cases, PCR duplicates were removed using Picard's MarkDuplicate (https://broadinstitute.github.io/picard/). Ancient DNA authentication was assessed using MapDamage [START_REF] Ginolhac | mapDamage: testing for damage patterns in ancient DNA sequences[END_REF] and contamination was estimated using both the mitochondrial genome and X chromosome [START_REF] Renaud | Schmutzi: estimation of contamination and endogenous mitochondrial consensus calling for ancient DNA[END_REF] .

Uniparental haplogroups

Consensus sequences of the mitochondrial genome were obtained using the software ANGSD v0.910, relying only on reads with base quality above 20, mapping quality above 30, and with more than 3X coverage. Mitochondrial haplotypes were determined based on the PhyloTree phylogenetic tree 42 build 17 and by using the HAPLOFIND web application [START_REF] Vianello | HAPLOFIND: a new method for high-throughput mtDNA haplogroup assignment[END_REF] and Phy-Mer [START_REF] Navarro-Gomez | Phy-Mer: a novel alignment-free and reference-independent mitochondrial haplogroup classifier[END_REF] . Unexpected or missing mutations were visually inspected in Geneious version R6 [START_REF] Kearse | Geneious Basic: An integrated and extendable desktop software platform for the organization and analysis of sequence data[END_REF] to check if they could be the results of misincorporations in low coverage regions. Sequences with more than 5% undetermined sites were removed. Y chromosome haplogroups were retrieved using the Yleaf software [START_REF] Ralf | Software for Human Y-Chromosomal Haplogroup Inference from Next-Generation Sequencing Data[END_REF] .

Nuclear markers

Target regions were summarized in a file, allowing their extraction using the mpileup tool from the SAMTools toolkit [START_REF] Li | A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data[END_REF] . Proper variant calling files were obtained using bcftools call and filtered for base quality 20 and a minimum coverage of 5 reads using the SelectVariants tool from the Genome Analysis Toolkit [START_REF] Mckenna | The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data[END_REF] . Allelic count and frequency estimation was carried out with plink [START_REF] Purcell | PLINK: a tool set for whole-genome association and population-based linkage analyses[END_REF] .

Reference datasets

We merged our dataset with published ancient individuals from Western Eurasia as well as genotyping data from the modern populations of the Human Origins panel [START_REF] Lazaridis | Genomic insights into the origin of farming in the ancient Near East[END_REF] . We only used published data obtained from UDG treated libraries. Genotype was called for the 1.2 million SNPs broadly used in aDNA in-solution capture procedures [START_REF] Mathieson | Genome-wide patterns of selection in 230 ancient Eurasians[END_REF] . A single allele was drawn at random for each position (minimum mapping and base quality of 30) using PileupCaller (https://github.com/stschiff/sequenceTools), rendering the individuals from the dataset homozygous for each locus. For genotyping data from modern individuals of the Human Origins panel, heterozygous sites in the eigenstrat file were randomly recoded as homozygous for the reference or the alternative allele using an awk script. The two files were then merged using plink v1.9 [START_REF] Purcell | PLINK: a tool set for whole-genome association and population-based linkage analyses[END_REF] .

  G2a and H2 have been described in various European Early Neolithic contexts in both Central Europe and Iberia[START_REF] Fu | The genetic history of Ice Age Europe[END_REF][START_REF] Olalde | Derived immune and ancestral pigmentation alleles in a 7,000-year-old Mesolithic European[END_REF][START_REF] Jones | Upper Palaeolithic genomes reveal deep roots of modern Eurasians[END_REF] , haplogroups I, whereas I2 are likely to have been introduced into the Neolithic pool through admixture with hunter-gatherers.In eastern France, the transition into the Middle Neolithic is represented by individuals from the Grossgartach culture, that derives from the Danubian sphere (~4,700-4,500 BCE, FigureS1) Their proximity to Central European contemporaneous individuals on PCA plots generated from both mitochondrial haplogroup frequencies and genome-wide markers suggests genetic homogeneity within the Danubian current (Figure1B, SI Appendix FigureS3-3). Interestingly, the BUCH2 individual harbors mixed GoyetQ2 and Villabruna HG ancestry (Figure1C) suggesting again that the GoyetQ2 ancestry is not a specific signature of south-western Europe. However, we can observe a subtle shift away from LBK individuals and towards WHG, reflecting an increase in shared alleles with the latter. This echoes observations from other ancient populations across Europe (Figure1C and 2A), yet with a higher proportion. This phenomenon seems to appear in France very early and becomes stronger as the Neolithic progresses, since individuals from French Middle Neolithic sites share a higher proportion of ancestry with HGs (Figure2A) and display Mesolithic signature haplogroups (e.g. U5b) at variable frequencies across the territory (see SI Appendix, FigureS3-2). Predominant in the Champagne region of northeastern France, these haplogroups suggest variable proportions of admixture between autochthonous hunter-gatherers and farmers. During the second half of the Middle Neolithic (MN2, ~4500 BCE), individuals from northern and southern France are poorly distinguished on the PCA plot (Figure 1B). Unlike in Iberia, the migration wave from which the Neolithic ancestry in ancient populations from France originated could not clearly be established using D statistics comparing the ancestry proportion between Iberia_EN and LBK_EN, as most individuals display non significant values (see SI Appendix, Figure S4-6) 12 . Using qpAdm, we investigated the ancestral source of the HG ancestry in individuals from French MN2 sites. Only individuals from southern France harbor varying levels of

  [START_REF] Olalde | The Beaker phenomenon and the genomic transformation of northwest Europe[END_REF] . French Beaker-associated individuals display a wide range of Steppe ancestry proportions (Figure1C, 2D).CBV95 in northern France derives the highest proportion of alleles from the Yamnaya in our dataset, and belongs to Y chromosome haplogroup R1b (Figure1C, 2D, SI Appendix FigureS4-5), providing the earliest clear evidence of the presence of this haplogroup in France around 2500 BCE (see SI Appendix, TableS5-1). This lineage was associated with the arrival of migrants from the Steppe in Central Europe during the Late Neolithic, and was described in other parts of Europe and in Bell Beaker associated individuals from southern France, while being almost absent in Iberia prior to the Bronze Age[START_REF] Olalde | The Beaker phenomenon and the genomic transformation of northwest Europe[END_REF][START_REF] Haak | Massive migration from the steppe was a source for Indo-European languages in Europe[END_REF] . PEI2, a male unearthed from a collective burial site near Carcassonne in southwestern France with artifacts of the Bell Beaker complex, falls within the genetic diversity of Neolithic individuals on the PCA. Modeling admixture proportions between three source populations:

Figure 1 :

 1 Figure 1: Overview of the ancient French dataset Filled symbols are used for ancient individuals from France whereas open symbols are used for other ancient Western Eurasian. The shape of the symbol indicates the geographic origin and the color, the time period. A. Location of the samples included in the study. B. Principal component analysis of ancient Western Eurasians projected onto the variation of present-day genotypes, restricted to Europe. C. Ancestry proportion for new and previously reported French individuals ranging from the Mesolithic to the Iron Age established using qpAdm. Each bar represents one individual with associated mtDNA haplogroup and Y chromosome haplogroup on the right side (error bars indicated S ±1SE). Symbol used in Figure 1A and B are indicated for each individual.

Figure 2 :C.

 2 Figure 2: Genetic affinity between ancient French and either Western hunter-gatherers or Steppe herders A. Boxplot of the D-statistics values of the form D(Mbuti, Test: WHG, Anatolia_Neolithic) for Early and Middle Neolithic European population with WHG consisting in Loschbour, ElMiron, GoyetQ2, Villabruna, KO1 and LaBrana to represent better the GoyetQ2 and Villabruna ancestry.

Figure 3 :

 3 Figure 3: Frequency of the derived allele for several types of genetic markers in Neolithic French and present-day Europeans. Numbers indicated above the bar chart correspond to the total number of observations of the derived allele for each variant in the Neolithic population. The color code reflects the function associated with the considered locus. The frequency of the alleles in present-day European populations is indicated in grey.

  I, although at a low frequency. First described in Eastern European populations and in Late Neolithic Central Europe[START_REF] Haak | Massive migration from the steppe was a source for Indo-European languages in Europe[END_REF] , these haplogroups suggest gene flow with a population deriving maternal ancestry from the Pontic Steppe herders at the onset of the Bronze Age. As foretold by its early occurrence in Bronze and Iron Age France share a common space on the PCA plot, both shifted towards modern Central Europe and falling within the genetic diversity of Bronze Age Britain and Central Europe (Figure1B), with a homogenization of the Steppe Component (Figure1C and 2C). In contrast to what was described for Central Europe[START_REF] Gamba | Genome flux and stasis in a five millennium transect of European prehistory[END_REF] , there is no further shift toward Eastern Eurasian genotypes during the Iron Age[START_REF] Gamba | Genome flux and stasis in a five millennium transect of European prehistory[END_REF] . Instead, the Steppe component, heterogeneously distributed between individuals during the Bronze Age (ranging from 30% to 70%), becomes homogeneous (Figure1Cand 2C), and individuals from the Hallstatt and La Tène culture on the French territory display similar affinities toward both modern and ancient populations. This could indicate that the transition from the Bronze Age to the Iron Age in France was mostly driven by cultural diffusion, without major gene-flow from an external population. This would be consistent with an archeological and linguistic hypothesis proposing that the Celts from the second Iron Age descended from populations already established in Western Europe, within the boundaries of the Bell Beaker cultural complex[START_REF] Brun | L'origine des Celtes. Communautés linguistiques et réseaux sociaux. The origin of the Celts. Linguistic communities and social networks[END_REF][START_REF] Demoule | Les indo-européens, un mythe sur mesure. Indo-europeans, a tailor-made myth La Recherche[END_REF][START_REF] Demoule | Mais où sont passés les Indo-Européens?: le mythe d'origine de l'Occident. But where have the Indo--europeans gone: the origin myth of the West[END_REF][START_REF] Hawkes | Celts, and Romans: studies in venture and resistance[END_REF] . It is important to mention, however, that due to the relative genetic homogeneity among European populations by the Bronze Age, subsequent migrations between different parts of Europe could easily remain unnoticed at this level of coverage

Beaker-associated CBV95, a drastic Y chromosome turnover occurs during the Bronze Age, where R1b replaces the preexisting diversity of Neolithic lineages in our sampling (Bronze Age: 11/13; Iron Age: 7/10 individuals having R1b, Table

S5

-1 and S5-2 and SI Appendix Figure

S5-1

).

Acknowledgments:

We would like to thank the Service Régional de l'Archéologie (SRA) from Occitanie, Grand Est, Hauts-de-France, ANTEA and INRAP for providing bone materials. We are grateful to Olivier Gorgé for making shotgun sequencing possible at the Institut de Recherche Biomédicale des Armées (IRBA) in the Division defense NRBC, département des services, Unité de Biologie moléculaire. Contributions: S.B., T.G., E.-M.G. and M.P. performed or supervised laboratory work with input from E.A.B., L.C., D.G.; archaeological materials were collected and information about archaeological contexts were provided by H.B.E., A.B., B.B., F.Ch., E.C., F.Co., B.D., S. D.-D., S.D., J.D., H.D., G.E., V.F., E.G., M.G., Y.G., S.G., J.G., L.H., M.I., F.L., F.M., B.P., S.P., E.P., I.P., I.R., V.R., R.R., B.S., C.T., S.T., E.V. and L.V.; S.B. and M.P analyzed data with input from E.A.B., T.G., E.-M.G.; S.B. and M.P. wrote the manuscript with the contribution of E.A.B., T.G., E.-M.G., J.D., B.B., M.G., J.G., R.R. and all other co-authors.

Funding: This work was supported by the Agence National de la Recherche (grant: ANR15-CE27-0001), the Institut National de Recherche en Archéologie preventive (INRAP) (PAS ANCESTRA) and the CNRS. The ARTEMIS program (French Ministry of Culture) funded the C14 datations. The paleogenomic facility of the Institut Jacques Monod (IJM) obtained support from the University Paris Diderot within the program "Actions de recherches structurantes". The sequencing facility of the IJM is supported by grants from the University Paris Diderot, the Fondation pour la Recherche Médicale (DGE20111123014), and the Région Ile-de-France (11015901).

Population Genetic Analyses

We conducted PCA on ancient individuals, projecting them onto the genetic variation of present-day Eurasians using the lsqproject option from smartpca [START_REF] Price | Principal components analysis corrects for stratification in genome-wide association studies[END_REF] . We computed D and f statistics using Admixtools to estimate gene flow and shared drift between populations and to test tree topologies [START_REF] Patterson | Ancient Admixture in Human History[END_REF] . We used qpAdm to estimate the ancestry proportions attributable to early Anatolian farmers, hunter-gatherers, and steppe herders in the individuals of this study [START_REF] Haak | Massive migration from the steppe was a source for Indo-European languages in Europe[END_REF] .

Data Availability:

Sequencing data can be found on the European Nucleotide Archive, ENA: ERP119731


Djeghri Nicolas
email: nicolas.djeghri@univ-brest.fr

Stibor Herwig

Lebeau Oanez

Pondaven Philippe

Nicolas Djeghri

δ13C, δ15N, and C:N ratios as nutrition indicators of zooxanthellate jellyfishes: insights from an experimental approach

Keywords: Stable isotopes, C:N ratios, Photosymbiosis, Mixotrophy, Scyphozoa, Zooxanthellae

Some jellyfish host zooxanthellae in their tissues (mostly from the family Symbiodiniaceae; Dinophyceae) and supplement their heterotrophic nutrition with their symbiont's photosynthates. The mixotrophy of zooxanthellate jellyfishes (as holobionts) renders the study of their nutrition, growth, and population dynamics complicated. Here, we used an experimental approach to assess how carbon and nitrogen stable isotopes (δ13C and δ15N) as well as the elemental composition (C:N ratios) of zooxanthellate jellyfishes are affected by variations in nutrition sources: i. e. predation (heterotrophic) versus photosynthesis (autotrophic). Our laboratory experiment, conducted on the zooxanthellate jellyfish Cassiopea sp. medusae (including symbionts) in the presence or absence of light and prey during 24 days, showed conclusive results. Presence of light decreased δ15N, increased δ13C and C:N ratios, whereas presence of prey increased δ15N, and decreased δ13C and C:N ratios. The medusae incubated with both light and prey had intermediate δ15N, δ13C and C:N ratios. Variations in zooxanthellate jellyfishes' nutrition sources (autotrophy vs. heterotrophy) are thus reflected by their isotopic and elemental composition. By disentangling the effects of autotrophy and of heterotrophy on zooxanthellate jellyfish isotopic and elemental compositions, these results would help to interpret the values of δ13C, δ15N and C:N ratios that can be observed on these organisms in fieldwork studies.

Introduction

Jellyfishes are increasingly acknowledged as an important component of marine ecosystems. Population dynamics of the pelagic life stages are often characterized by important fluctuations with dramatic biomass increases followed by sudden collapses (Lucas andDawson 2014, Pitt et al. 2014). These fluctuations can have important consequences for pelagic community dynamics and nutrient cycling (Pitt et al. 2009a), or for human activities [START_REF] Purcell | Anthropogenic causes of jellyfish blooms and direct consequences for humans: a review[END_REF]). One of the key factors controlling jellyfish population dynamics, is nutrition (e. g. Lucas and Dawson 2014, [START_REF] Pitt | Bloom and bust: Why do blooms of jellyfish collapse?[END_REF]). One way to study jellyfish nutrition is to use their stable isotopes signatures (mainly δ 13 C and δ 15 N, see Pitt et al. 2009b).

Many recent studies have focused on jellyfish stable isotopes and have provided precisions of their diets as well as competition relationships (e. g. Fleming Zooxanthellate jellyfishes (mostly Rhizostomeae, Scyphozoa, see [START_REF] Djeghri | Review of the diversity, traits, and ecology of zooxanthellate jellyfishes[END_REF]) are characterized by their photosymbiotic relationship with zooxanthellae (generally from the family Symbiodiniaceae, Dinophyceae; LaJeunesse 2001, [START_REF] Lajeunesse | Systematic revision of Symbiodiniaceae highlights the antiquity and diversity of coral endosymbionts[END_REF]). This symbiotic relationship is thought to be similar to the one well known in corals with the zooxanthellae providing their host with photosynthates while recycling the host's respiration and excretion products (see [START_REF] Davy | Cell biology of Cnidarian-Dinoflagellate symbiosis[END_REF]). In such symbiosis, both carbon and nitrogen can be obtained via heterotrophy and autotrophy and are recycled between the host and its zooxanthellae. Carbon dioxide (CO2) and nitrogen from host respiration and excretion are used and metabolized by the zooxanthellae. Simultaneously, complex molecules (including carbohydrates, lipids and amino acids) are transferred from one partner to the other [START_REF] Davy | Cell biology of Cnidarian-Dinoflagellate symbiosis[END_REF]).

Zooxanthellate jellyfishes, as holobionts (host + symbionts), are thus mixotrophs, deriving their nutrition from both predation and zooxanthellae's photosynthesis (Kremer 2005, Welsh et al. 2009). Generally, the symbiosis provides most if not all of the carbon needed for respiration [START_REF] Kremer | Significance of photosynthetic endosymbionts to the carbon budget of the schyphomedusa Linuche unguiculata[END_REF][START_REF] Kikinger | Cotylorhiza tuberculata (Cnidaria: Scyphozoa) -Life history of a stationary population[END_REF], McCloskey et al. 1994, Verde and McCloskey 1998) while predation is still needed to meet nitrogen and phosphorus requirements (Kremer 2005, Welsh et al. 2009). However, the relative contribution to nutrition of the predation versus the photosynthesis might be variable across species, populations, environments, or during growth (see e. g. [START_REF] Zeman | Trophic ecology of the neustonic cnidarian Velella velella in the northern California Current during an extensive bloom year: insight from gut contents and stable isotope analysis[END_REF]). The conclusions of these studies have, however, been limited due to the lack of data on the interplay between autotrophy and heterotrophy of zooxanthellate jellyfishes as reflected in their isotopic and elemental composition [START_REF] Zeman | Trophic ecology of the neustonic cnidarian Velella velella in the northern California Current during an extensive bloom year: insight from gut contents and stable isotope analysis[END_REF]). To better understand this, controlled experiments are needed where the resources for heterotrophy (prey), and for autotrophy (light) can be manipulated and their effect on stable isotopes signatures and elemental composition can be assessed. In this study, we aim to provide experimental insights on how isotopic and elemental composition of zooxanthellate jellyfishes are affected by different food regimes, and, more specifically, by variations of the relative importance of autotrophy and heterotrophy. In order to achieve this, we assessed the changes in the δ 13 C, δ 15 N, and C:N ratios in young specimens of zooxanthellate Cassiopea sp. medusae (Scyphozoa: Rhizostomeae) over a period of 24 days and in the presence or absence of prey and light.

Materials and Methods

Experimental set-up

Small specimens of Cassiopea sp. medusae (ca. 6 mm in bell diameter and 1 month old) were acquired from the Trocadéro Aquarium (Paris, France). In this aquarium, the medusae were kept at 25 °C with a daily light cycle, and were fed Artemia sp. nauplii twice a day. After their arrival to the laboratory the medusae were acclimatized to local heated (25-26 °C) filtered (1 μm pore size) sea-water during one day.

In the following day, five medusae were randomly sampled to represent the initial state and then, the experiment was set up. A total of 72 medusae specimens were individually placed in 75 ml glass flasks filled with 50 mL of filtered sea-water. The flasks were then randomly assigned to one of the four experimental treatments (18 medusae per treatment). The experimental treatments were as follows: (1) fed and in light, (2) fed and in the dark, (3) starved and in light, and (4) starved and in the dark. The goal of these different treatments was to target respectively: mixotrophy, heterotrophy, and autotrophy, the fourth treatment being a control. The flasks containing the medusae were kept in water baths, which regulated a fixed temperature (25-26 °C). Two water baths were used, one for the medusae kept in light, and the other for the medusae kept in the dark. The temperature changed little during the experiment and between the two water baths (25.6 ± 0.4 °C and 25.3 ± 0.4 °C respectively in the lighted and darkened water baths; mean ± s.d.). The light was provided by a fluorescent lamp on a 12:12 hours day:night cycle at an intensity of ca. 110 µmol photons.m -2 .s -1 . Food consisted of 2 h of feeding ad libitum every two days on Artemia sp. young nauplii (< 24 h after hatching). The medusae full guts and numerous remaining prey in the flasks confirmed a proper ad libitum feeding. Every two days, and after the feeding, the incubation water of the medusae was entirely changed. During the latter process, the flasks were also washed to prevent fouling.

Every four days, and before the feeding (ensuring empty guts), three medusae were randomly sampled from each treatment. Thus the sampling occurred at the following days: 4, 8, 12, 16, 20 and 24; plus the initial state being represented by the five medusae sampled before setting the treatments.

Processing of medusae

Immediately after sampling, presence and physiological state of the symbionts were assessed. For this purpose, the medusae were put in the dark for at least 20 minutes allowing the opening of the photosystem reactive centers of zooxanthellae. The photosynthetic parameters of the zooxanthellae were assessed by pulse amplitude modulatory (PAM) fluorometry using the "OJIP protocol" of an AquaPen-C-AP-C100 (®Photon Systems Instruments, PSI, Brno, Czech Republic), at a 450 nm excitation wavelength. This protocol measures the fluorescence emitted after a flash of saturating light. This yields the estimation of several variables among which the maximum photosystem II quantum yield (Fv/Fm), which is a proxy of the photosystem II efficiency. With Fm the maximum fluorescence under saturating light, and Fv = Fm -F0 with F0 the initial fluorescence (Strasser et al. 2000). The measures were performed on whole medusae specimens. Prior to the measurement, it was ensured that the medusae were settled at the cuvette bottom to insure proper exposition to the saturating flash. Additionally, along each medusae measure, a blank was realized by using the "OJIP protocol" on incubation water without medusae (75 blanks total).

Preparation for elemental and isotopic analysis started first by quickly rinsing the medusae in deionized water in order to remove the sea salt. The whole medusae were then placed in pre-weighted tin capsules (10.5 × 9 mm EMAL technology, United Kingdom) and oven-dried at 60 °C for ca. 48 h. After drying, the tin capsules containing the medusae were locked and weighted again to obtain the medusae dry mass, which varied between 0.2 and 3.2 mg (0.9 ± 0.6 mg; mean ± s.d.). In the preparation procedure, the medusae were unpreserved before the isotopic analysis, following the recommendations of [START_REF] Fleming | Preservation methods alter stable isotope values in gelatinous zooplankton: implications for interpreting trophic ecology[END_REF]. Due to the small size of the medusae specimens, it was not possible in this experiment to follow the recommendations of MacKenzie et al. ( 2017) by dissecting, washing and rubbing the mesoglea before conservation. In addition, due to their small size, it was not possible to separate the animal tissue from the zooxanthellae. This would have resulted to acquiring an insufficient biomass for the isotopic analysis. Thus, the results presented here are measures on the holobiont (animal host + zooxanthellae).

Processing of prey nauplii

In order to assess the isotopic and elemental composition of the Artemia sp. nauplii given as food for the medusae, we sampled them two times during the experiment (at days 8 and 14) plus a sampling at day 28, slightly after the end of the experiment (protocol unchanged). At each sampling, concentrated nauplii were divided in five aliquots, and oven-dried at 60 °C for ca. 48 h in clean glass flasks. The dried nauplii were then scratched from the flasks and ground into a powder. Finally, between 0.5 and 1.5 mg of the powder were inserted and locked in tin capsules (10.5 × 9 mm EMAL technology, United Kingdom).

Elemental and stable isotopes composition

The analyses of medusae and nauplii samples were performed using an Elemental analyzer (Thermo Scientific EA Flash 2000), coupled to a Mass Spectrometer (Thermo Scientific DELTA V Plus) at the Stable Isotopes Laboratory of the "Pôle Spectrométrie Océan" (PSO-IUEM, Plouzané, France). The nitrogen and carbon mass of medusae samples ranged respectively from 15 to 109 µgN (35 ± 22 µgN; mean ± s.d.) and from 60 to 543 µgC (175 ± 120 µgC; mean ± s.d.). As the whole medusae were inserted in the tin capsules, these values are representative of their total weights. The nitrogen and carbon mass of nauplii samples ranged respectively from 40 to 105 µgN (64 ± 18 µgN; mean ± s.d.) and from 200 to 522 µgC (326 ± 90 µgC; mean ± s.d.). The samples were calibrated for mass bias using casein (IVA-33802155, Analysentechnik, Germany) as the elemental standard (range: 5-108 µgN; 16-377 µgC). Some materialrich samples were automatically diluted during the analysis process (Thermo Scientific ConFlo IV).

Stable isotopes values are expressed as permil (‰) using the δ notation (normalized to Vienna Pee Dee Belemnite and atmospheric N2 for respectively carbon and nitrogen):

δX = ( X sample H X sample L ⁄ X std H X std L ⁄ -1) × 1000
With X the element measured, X H the amount of the heavy isotope and X L the amount of the light isotope from the samples (Xsample) and the standard (Xstd).

As some of our samples had a low (< 20 µgN) nitrogen mass, we analyzed five replicates of casein standards with a low nitrogen mass (13.4 ± 1.9 µgN; mean ± s.d.) to check whether this low mass have led to uncertainties in our measures. We found only a low variability on the obtained δ 15 N measures (0.04 ‰ s.d., n=5) indicating that our measures were consistent even at low biomass levels.

Unless indicated otherwise, all C:N ratios are expressed by mass (following Ikeda 2014 and Molina-

Ramírez et al. 2015)

. As the C:N ratios of both the medusae and their prey were higher than 3.5, a normalization of the δ 13 C for lipid content was advisable [START_REF] Post | Getting to the fat of the matter: models, methods and assumptions for dealing with lipids in stable isotopes analyses[END_REF]). For the nauplii, we used the general normalization for aquatic animals given by [START_REF] Post | Getting to the fat of the matter: models, methods and assumptions for dealing with lipids in stable isotopes analyses[END_REF] and for the medusae, we used the normalization specific to scyphozoans proposed by D'Ambra et al. (2014). For comparison, raw data is still presented as supplementary material (see discussion).

Statistics

The data collected during the experiment (carbon masses, δ 13 C, δ 15 N, C:N ratios and Fv/Fm) were analyzed using linear mixed-effects models (LME) (e. g. [START_REF] Crawley | The R book[END_REF]). The presence or absence of prey and light were considered as fixed effects while time was considered as random effect. Model assumptions (mean of residuals = 0, linearity and normality) were checked using model-checking plots. If the model assumptions were not met, the data were log transformed. If the fixed effects affected significantly the results (if p-value < 0.05), subsequent Tukey post-hoc tests were performed on leastsquare means (i. e. means adjusted for the effect of time) to assess which combination of the fixed effects (light and prey) led to different responses.

One-way ANOVAs were used to assess possible variations in prey δ 13 C, δ 15 N, and C:N ratios over time.

Normality and homoscedasticity assumptions were checked using the Shapiro-Wilk normality test and Bartlett homogeneity of variance test (threshold: α = 0.05). If the assumptions were not met, the data were Box-Cox transformed. All statistical analyses were performed in R (R Core Team 2017).

Results

Mass variation in medusae

The carbon mass of the medusae (Fig. 3.2. δ 13 C and δ 15 N

The medusae δ 13 C values (Fig. 2 A,B) were significantly affected by light (LME, t-value = 22.5, p-value < 0.001) and by the interaction of light and prey (LME, t-value = -6.2, p-value < 0.001). At the beginning of the experiment δ 13 C values of the medusae were of -18.7 ± 0.9 ‰ (mean ± s.d.). In the treatment with both light and prey the δ 13 C values of the medusae increased quickly (in less than four days) reaching ca.

-15 ‰. This trend was even more pronounced in the treatment with only light were the δ 13 C values of the medusae reached ca. -13 ‰. Conversely, in the treatments with only prey or without resources the δ 13 C values of the medusae remained similar or decreased slightly throughout the experiment (generally comprised between -19 ‰ and -21 ‰). It should be noted that the distinction between the medusae from the treatments with light alone and with light and prey is not as distinct with data not normalized for lipids (Fig. S1).

The medusae δ 15 N values (Fig. 2

C:N ratios

The C:N ratios obtained in the medusae (Fig. 3) were significantly affected by light (LME, t-value = 27.5, p-value < 0.001) and by the interaction of light and prey (LME, t-value = -6.6, p-value < 0.001). At the beginning of the experiment, the C:N ratios of medusae were of 4.5 ± 0.1. These values increased in the treatment with both light and prey (4.8 ± 0.1 at the end of experiment; mean ± s.d.). A similar, but more pronounced increased was seen in the C:N ratios of the medusae exposed to only light (5.9 ± 0.2 at the end of experiment; mean ± s.d.). By opposition, C:N ratios decreased slightly in the treatments with only prey or without resources (respectively reaching 4.1 ± 0.1 and 3.9 ± 0.1 at the end of the experiment; mean ± s.d.). The C:N ratios obtained in Artemia sp. nauplii did vary significantly during the experiment (ANOVA, F = 25.9, p-value < 0.001). However, these variations were of small amplitudes (minimum: 4.97, maximum:

5.20) compared to the variationsobserved in Cassiopea sp. medusae following the different experimental treatments (Fig. 3 A). Therefore, these small variations of the C:N ratios of the prey are unlikely to have significantly affected the outcome of the experiment. Throughout the experiment, C:N ratios in Artemia sp. nauplii averaged 5.1 ± 0.1 (mean ± s.d.).

PAM parameters of zooxanthellae

The blanks always yielded low values of F0 (90 ± 9; mean ± s.e.m.) as compared to the F0 values of the medusae (5270 ± 630; mean ± s.e.m.). This equates to a signal-to-noise ratio of ca. 60, which is sufficient to have a reliable estimate of photosynthetic activity. Two outliers were removed from the medusae's PAM data (Fv/Fm below 0.4, similar to a blank, most likely due to a lack of exposition of the medusae to the saturating flash). With the exception of this two outliers, the Fv/Fm of medusae remained very stable in all conditions and during the whole experiment averaging an overall value of 0.70 ± 0.06 (mean ± s.d.;

Fig. S2). The LME models did not indicate any effect of presence or absence of prey and light on the zooxanthellae Fv/Fm. Independently of this lack of effect of the experimental treatments on PAM parameters (discussed in Supplementary Material 2), high F0 values as compared to the blanks demonstrate the presence of zooxanthellae in the Cassiopea sp. medusae in all treatments and during the whole experiment.

Discussion

Isotopic composition

In this study, the δ 13 C values obtained in the medusae were the highest in the treatment with light only, lowest in the treatments with only prey and without resources, and intermediate in the treatment with both prey and light (Fig. 2A,B). Similar effects of heterotrophic feeding on δ 13 C values have been reported for corals (e. g. [START_REF] Reynaud | Effect of feeding on the carbon and oxygen isotopic composition in the tissues and skeleton of the zooxanthellate coral Stylophora pistillata[END_REF][START_REF] Ferrier-Pagès | Summer autotrophy and winter heterotrophy in the temperate symbiotic coral Cladocora caespitosa[END_REF]). As in this study, the δ 13 C of the predator tended towards the δ 13 C of the prey when fed. However, some caution should be taken when interpreting the results of this study concerning the effects of heterotrophic feeding on zooxanthellate jellyfish's δ 13 C. Indeed, it is unsure that the lipid normalization used here can be applied to a photosymbiotic holobiont as it has been derived from the heterotrophic Aurelia sp. (D'Ambra et al.

2014

). Without this normalization, the effect of heterotrophic feeding on δ 13 C is less clear (Fig. S1). Thus, albeit an effect of heterotrophic feeding on zooxanthellate jellyfishes' δ 13 C is likely, our results should be taken with caution regarding this point. In contrast, light had a clear positive effect on δ 13 C of zooxanthellate jellyfishes whether the normalization for lipid content is made or not (Fig. 2A,B; Fig. S1).

These conclusions are consistent with the previous findings on corals (e. g. In the experiment, δ 15 N values were the lowest in the medusae exposed to only light, compared to the other treatments (Fig. 2C,D). This is different from what is known in tropical scleractinian corals in which photosynthesis tend to increase, or have little effect on δ 15 N rather than decrease it, like seen here (see with no resources and the treatment with only prey presented the same δ 15 N (Fig. 2C,D). However, the effect of predation is clear as the δ 15 N values in the treatment with prey and light were intermediate between those of the treatment with only light, and the treatment with only prey (Fig. 2C,D). The similitude between the treatment without resources and the treatment with only prey would thus be explained by the initial condition (i. e. at day 0, medusae already had high δ 15 N). Thus, overall, predation would have led to higher δ 15 N of medusae (Fig. 2C,D). Interestingly however, the δ 15 N values in medusae were never higher than those obtained on prey (Fig. 2C), which suggests that no measurable fractionation occurred between the holobiont and their prey. This is most likely due to high recycling of nitrogen between the host and its symbionts (see also [START_REF] Reynaud | Effect of light and feeding on the nitrogen isotopic composition of a zooxanthellate coral: role of nitrogen recycling[END_REF]).

These patterns seen in δ 13 C and δ 15 N can be explained through two main processes affecting stable isotopes signatures in photosymbiotic organisms (reviewed in Ferrier-Pagès and Leal 2018):

The first process is the mixing of carbon or nitrogen coming from two contrasted sources; autotrophic uptake of dissolved inorganic nutrients, on one hand, and heterotrophic predation-mainly Thus, δ 13 C values obtained through zooxanthellae's autotrophy would be higher than those obtained through predation on zooplankton (Fig. 4A). For nitrogen, the pattern is reversed; zooxanthellae take up dissolved inorganic nitrogen with a low δ 15 4B). For both carbon and nitrogen, the isotopic signature of the two sources (inorganic nutrient uptake, and predation) is then exchanged and recycled between the zooxanthellae and the host (e. g. [START_REF] Reynaud | Effect of light and feeding on the nitrogen isotopic composition of a zooxanthellate coral: role of nitrogen recycling[END_REF]).

The second process involves the depletion of in-hospite nutrient pools due to photosynthesis.

Zooxanthellae tend to take up preferentially inorganic nutrients with light isotopes resulting in fractionation (Ferrier-Pagès and Leal 2018). However, at high photosynthesis rates, the host's pool of inorganic nutrients can get depleted. Thus, to meet their photosynthetic requirements, zooxanthellae take up more heavy isotopes, reducing fractionation ("depletion-diffusion hypothesis", see Muscatine et al. 1989a, Fig. 4A andB). This results in a tendency for isotopic signature to correlate with photosynthesis levels. The higher the photosynthesis, the higher the δ 13 It is important to notice that these two processes-mixing of the heterotrophic and autotrophic sources, and reduced fractionation at high photosynthesis levels-would have similar consequences on δ 13 C, but not on δ 15 N. For δ 13 C, a predominantly autotrophic nutrition would imply a decreased fractionation inhospite due to high photosynthesis and a predominant uptake of dissolved inorganic nutrients. Both of these processes result in high δ 13 C (and vice-versa for a predominantly heterotrophic nutrition, Fig. 4A).

By opposition, for δ 15 N, reduced fractionation due to high photosynthesis levels would result in high δ 15 N, but high uptake of dissolved inorganic nutrients would result in low δ 15 N (and vice-versa for a predominantly heterotrophic nutrition, Fig. 4B). Thus, to understand how δ 15 N would react to change in holobiont nutrition, it is important to know which of the above-mentioned processes controls its dynamics.

In the context of this study, the increase in δ 13 C due to light, could be explained by both increased uptake of CO2 and reduced fractionation at higher photosynthesis (and conversely in the dark, Fig. 4A).

The results obtain on δ 15 N values in this study, suggest that, in zooxanthellate jellyfishes, of the two processes above-mentioned-mixing of autotrophic and heterotrophic sources, and reduced fractionation at high photosynthesis-the former is the dominant one (Fig. 4B). Hence, in zooxanthellate jellyfishes, a predominantly autotrophic nutrition would imply that most nitrogen comes from the fixation of dissolved inorganic nitrogen (see e.

C:N ratios

In our experiment, light increased C:N ratios of the Cassiopea sp. holobionts whereas prey generally decreased it, with intermediate C:N ratios in the treatment with both light and prey (Fig. 3A,B). As with δ 15 N, C:N ratios of the treatment with no resources and with prey only were similar. This similarity may be again explained by the initial condition. Another potential explanation would be that starved zooxanthellate jellyfishes use first reserves accumulated by photosynthesis (generally carbon rich, e. g.

Muller-Parker et al. 1996). Such preferential degradation would make their C:N ratios decrease and get similar to the ones typically reported for non-zooxanthellate jellyfishes [START_REF] Ikeda | Synthesis toward a global model of metabolism and chemical composition of medusae and ctenophores[END_REF][START_REF] Molina-Ramírez | Functional differences in the allometry of the water, carbon and nitrogen content of gelatinous organisms[END_REF]. Independently of the treatment with no resources, our results suggest that predation would tend to decrease C:N ratios (Fig. 3A,B). Such a decrease of C:N ratios due to predation have already been reported for the zooxanthellae of a sea anemone [START_REF] Cook | Host feeding and nutrient sufficiency for zooxanthellae in the sea anemone Aiptasia pallida[END_REF]). However, other studies have pointed out that a similar decrease of C:N ratios can also be due to an enrichment by dissolved inorganic nitrogen (e. g. is likely that their C:N ratios would react to nutrient enrichment too. This suggests that C:N ratios of zooxanthellate jellyfishes might be impacted by nitrogen availability (either as prey or as dissolved inorganic nitrogen).

Remarks on tissue turnover

One of the advantages of the study of the elemental and isotopic composition over e.g. gut content analyses, is that it provides a more time-integrated information (Pitt et al. 2009b). The time frame represented by isotopic composition is, however, dependent on tissue turnover which can be variable as function of, e. g. taxonomy, organ, body size, or temperature [START_REF] Thomas | Predicting rates of isotopic turnover across the animal kingdom: a synthesis of existing data[END_REF]Crowther 2015, Vander Zanden et al. 2015). In another scyphozoan jellyfish, Aurelia sp., the isotopic half-life was determined to be ca. 10 days for both carbon and nitrogen (D'Ambra et al. 2014). In the present experiment, changes of δ 13 C and C:N ratios occurred very fast (within the first four days, Figs. 2A,3). This was apparently less true for the δ 15 N values which may have experienced slower changes (Fig. 2B). These fast changes may have several explanations: First, the medusae used here were of small size and, additionally, some of them grew (Fig. 1, see also Supplementary Material 2) which can explain the fast changes (

Implications for fieldwork studies

One of the challenges to understand the nutrition of zooxanthellate jellyfishes in their natural environments relates to their mixotrophy. As zooxanthellate jellyfishes obtain their nutrition from predation and photosynthesis (Kremer 2005, Welsh et al. 2009), both processes must be investigated.

Ideally, predation, photosynthesis, respiration, nutrient uptake and excretion have all to be measured which may represent an intensive amount of work rarely carried out in its entirety (see however, [START_REF] Kremer | Significance of photosynthetic endosymbionts to the carbon budget of the schyphomedusa Linuche unguiculata[END_REF], Kremer 2005). Studies of stable isotopes and elemental composition are comparatively easier and have the advantage of providing more time-integrated information (Pitt et al. 2009b). The findings of this study provide baseline information on how C:N ratios, δ 13 C and δ 15 N can be interpreted in fieldwork studies focusing on the nutrition sources of zooxanthellate jellyfishes.

To summarize, our results suggest that both δ 13 C and δ 15 N vary with the relative balance of autotrophy vs. heterotrophy (although, the effect of heterotrophy on δ 13 C is less well supported by our results). It would be expected that, if dominantly heterotrophic, zooxanthellate jellyfishes would have δ 15 N values close to those of their prey (see e. g. [START_REF] Zeman | Trophic ecology of the neustonic cnidarian Velella velella in the northern California Current during an extensive bloom year: insight from gut contents and stable isotope analysis[END_REF]. By opposition, if dominantly autotrophic, zooxanthellate jellyfishes would have δ 15 N values close (or lower) than those of primary producers (see e. g. [START_REF] Freeman | Variation in δ 13 C and δ 15 N values suggests a coupling of host and symbiont metabolism in the Symbiodinium-Cassiopea mutualism[END_REF]). The δ 13 C values would display opposite trends. Finally, C:N ratios may be indicators of the efficiency of nitrogen supplies. Future fieldwork studies would be able to build on these results to better characterize zooxanthellate jellyfishes' nutrition.

Funding

  ) was significantly affected by light alone (LME, t-value = 7.5, pvalue < 0.001) and by the interaction between light and prey (LME, t-value = 3.4, p-value < 0.01). At the beginning of the experiment, the carbon mass of Cassiopea sp. medusae was of 115.1 ± 27.0 µg C (mean ± s.d.). Only the carbon mass of the medusae in the treatment with both light and prey did noticeably increased, reaching 449.3 ± 68.2 µg C (mean ± s.d.) at the end of the experiment. The carbon mass of the medusae in the treatment with only light did not increased significantly reaching a carbon mass value of 133.9 ± 27.9 µg C (mean ± s.d.) at the end of the experiment. The carbon mass values of the medusae in the treatments with only prey or without resources tended to decrease, dropping to respectively 88.0 ± 26.1 µg C and 71.1 ± 15.7 µg C (mean ± s.d.) at the end of the experiment. In the treatment without resources, one medusae specimen died. Thus for this treatment, there remained two replicates instead of three at day 24.

Fig. 1 .

 1 Fig. 1. (A) Changes in the Cassiopea sp. medusae carbon mass (µg C) (means ± s.e.m.) over the course of the experiment as a function of the experimental conditions. (B) Comparison of least-square means obtained from each treatment (± 95 % C. I.). The letters (a, b, and c) indicate statistically different treatments (Tukey post hoc test, p-value < 0.05).

  C, D) were significantly affected by both light (LME, t-value = -8.7, pvalue < 0.001) and prey (LME, t-value = 2.2, p-value < 0.05), but not by their interaction. At the beginning of the experiment, the δ 15 N values of medusae was of 8.9 ± 1.1 ‰ (mean ± s.d.). These values decreased slightly in the treatment with light and prey reaching 8.0 ± 0.3 ‰ (mean ± s.d.) at the end of the experiment. The decrease was more pronounced in the treatment with light only which reached 5.5 ± 0.4 ‰ at the end of the experiment. Finally, δ 15 N values did not change in the treatments with only prey or without resources (values at the end of the experiment of 9.0 ± 0.3 ‰ and 8.8 ± 1.0 ‰ respectively; mean ± s.d.). The δ13 C and δ15 N obtained in Artemia sp. nauplii prey did not vary significantly during the experiment (ANOVAs, p-values > 0.05) averaging respectively -19.4 ± 0.2 ‰ and 10.7 ± 0.5 ‰ (mean ± s.d.) (Fig.2 A, C).

Fig. 2 .

 2 Fig. 2. Changes in the Cassiopea sp. medusae δ 13 C (‰) (A) and δ 15 N (‰) (C) (means ± s.e.m.) over the course of the experiment as a function of the experimental conditions. Solid and dashed blue lines represent the mean ± s.d. of the isotopic signatures of Artemia sp. nauplii used as prey in fed treatments. (B and D) Comparison of least-square means obtained for each treatment (± 95 % C. I.). The letters (a, b, and c) indicate statistically different treatments (Tukey post hoc test, p-value < 0.05). δ 13 C values of have been normalized for lipid content according to Post et al. (2007) for nauplii, and D'Ambra et al. (2014) for Cassiopea sp. medusae.

Fig. 3 .

 3 Fig. 3. (A) Changes in the Cassiopea sp. medusae mass C:N ratios (means ± s.e.m.) over the course of the experiment as a function of the experimental conditions. The blue error bar indicate the range of C:N ratios in nauplii prey measured in the course of the experiment (see text). (B) Comparison of the leastsquare means obtained for each treatment (± 95 % C. I.). The letters (a, b, and c) indicate statistically different treatments (Tukey post hoc test, p-value < 0.05).

  on zooplankton-on the other hand(Reynaud et al. 2002, Alamaru et al. 2009, Ferrier-Pagès et al. 2011, Ferrier-Pagès and Leal 2018). The uptake of dissolved inorganic carbon by zooxanthellae generally leads to higher δ 13 C values (typically -10 ‰ to -14 ‰) than those of typical oceanic particulate organic matter and plankton (ca. -20 ‰; Muscatine et al. 1989a, Ferrier-Pagès et al. 2011, Ferrier-Pagès and Leal 2018).

  N value (ca. 5 ‰ Ferrier-Pagès and Leal 2018) while predation leads the uptake of nitrogen with higher δ 15 N values due to fractionation through the food web (Post 2002, Ferrier-Pagès et al. 2011, Ferrier-Pagès and Leal 2018, Fig.

  photosynthesis levels. The higher the photosynthesis, the higher the δ13 C (Muscatine et al. 1989a, Swart et al. 2005, Alamaru et al. 2009) or the δ 15 N values (Muscatine and Kaplan 1994, Baker et al. 2011; reviewed in Ferrier-Pagès and Leal 2018).

  g.[START_REF] Muscatine | Dissolved inorganic nitrogen flux in symbiotic and nonsymbiotic medusae[END_REF][START_REF] Wilkerson | DIN, DON and PO4 flux by medusa with algal symbionts[END_REF][START_REF] Freeman | Symbiont carbon and nitrogen assimilation in the Cassiopea-Symbiodinium mutualism[END_REF], and would result in low δ15 N. On the contrary, a predominantly heterotrophic nutrition would imply that more nitrogen comes from predation (mainly on zooplankton) resulting in a comparatively higher δ 15 N. Thus, values of δ15 N can be considered as a good indicator of the relative importance of autotrophy and heterotrophy in zooxanthellate jellyfishes.

Fig. 4 .

 4 Fig. 4. Conceptual diagram illustrating how δ 13 C (A) and δ 15 N (B) are affected by the relative importance of heterotrophic (blue circles and arrows, Chetero and Nhetero) and autotrophic (orange circles and arrows, Cauto and Nauto) nutrition pathways in zooxanthellate jellyfishes. DIC = dissolved inorganic carbon, DIN = dissolved inorganic nitrogen, Phyto = phytoplankton, Zoo = zooplankton, Phot = photosynthesis, Frac = fractionation (more photosynthesis tend to decrease fractionation, this effect, in zooxanthellate

  [START_REF] Fry | Rapid 13 C/ 12 C turnover during growth of brown shrimp (Penaeus aztecus)[END_REF][START_REF] Thomas | Predicting rates of isotopic turnover across the animal kingdom: a synthesis of existing data[END_REF][START_REF] Vander Zanden | Stable isotopes turnover and half-life in animal tissues: a literature synthesis[END_REF]. In its natural environment, Cassiopea sp. can grow up to ca. 20-25 cm in bell diameter (see e. g.[START_REF] Morandini | All non-indigenous species were introduced recently? The case study of Cassiopea (Cnidaria: Scyphozoa) in Brazilian waters[END_REF]. It is unlikely that such large specimens would display such fast change in composition. Another aspect that could explain the fast change in elemental and isotopic composition observed here, is that Cassiopea sp. is zooxanthellate. The zooxanthellae are also likely to impact residence time of elements within the holobiont, possibly differently for nitrogen and carbon, due to recycling[START_REF] Reynaud | Effect of light and feeding on the nitrogen isotopic composition of a zooxanthellate coral: role of nitrogen recycling[END_REF]).

  et al. 2015, Javidpour et al. 2016, Vansteenbrugge et al. 2016, D'Ambra et al. 2018, Milisenda et al. 2018). Most of these studies have strictly focused on heterotrophic jellyfishes. However, some jellyfishes are known to live in symbiosis with zooxanthellae. Comparatively, the zooxanthellate jellyfishes have received less interest (see however Freeman et al. 2017, Zeman et al. 2018).

  [START_REF] Sugiura | On the life-history of Rhizostome medusae V. On the relation between zooxanthellae and the strobilation of Cephea cephea[END_REF][START_REF] Mccloskey | Daily photosynthesis, respiration, and carbon budgets in a tropical marine jellyfish (Mastigias sp.)[END_REF], Verde and McCloskey 1998, Bolton and Graham 2004, reviewed in Djeghri et al. 2019). Studies using stable isotopes, in this context might be valuable tools to understand these variations. Numerous studies on other photosymbiotic cnidarians (mainly corals), have shown that variations of nutrition affect the isotopic and elemental composition (see e. g. Muscatine et al. 1989a, Muscatine and Kaplan 1994, Alamaru et al. 2009, Reynaud et al. 2009, Ferrier-Pagès et al. 2011, reviewed by Ferrier-Pagès and Leal 2018). Similar effects can be expected in zooxanthellate jellyfishes. To date, only few fieldwork studies have focused on the isotopic composition of zooxanthellate jellyfishes (see Freeman et al. 2017,

  Muscatine et al. 1989a, Swart et al. 2005, Alamaru et al. 2009, Ferrier-Pagès et al. 2011).

  Muscatine and Kaplan 1994, Alamaru et al. 2009, Reynaud et al. 2009). Our results are more comparable to what is observed in more heterotrophic temperate corals (Ferrier-Pagès et al. 2011). The treatment

This work is part of the Ph.D. research program of ND, funded by the Université de Bretagne Occidentale (UBO, Brest, France). This work was also supported by the Labex Mer (ANR-10-LABX-19-01, IUEM, Brest, France).

Acknowledgements

The authors would like to thank Eric Dabas, Leslie Gager and William Handal for providing technical assistance during the experiments. We also thank Gauthier Schaal and Rudolph Corvaisier for their valuable comments on the first version of this manuscript. Finally, we wish to thank two anonymous reviewers for their helpful comments.

Supplementary Material 2: Notes on growth of medusae and PAM parameters

Medusae grew only in the treatment with both light and prey (i. e. increase in carbon mass of the holobiont, Fig. 1 in main text; and increase in medusae bell diameter from ca. 6 mm to 8-10 mm, data not shown). The absence of carbon mass increase (Fig. 1