Marta Gallardo Ruiz

Jean-François Le Galliard

Thomas Tully

Marta Gallardo Ruiz
email: 144274730.marta.gallardo_ruiz@etu.upmc.fr

Marta Gallardo

Ruiz +33

Genetic variation in light vision and light-dependent movement behaviour in the eyeless Collembola Folsomia candida

Keywords: Folsomia candida, Collembolan, light sensitivity, eyeless species, mobility, intraspecific variation

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

Introduction

Animals use various environmental cues for habitat choices, and different individuals from a single species may show contrasting habitat preferences depending on their sex, stage, size and genetic background for instance [START_REF] Stamps | Habitat selection by dispersers: integrating proximate and ultimate approaches[END_REF][START_REF] Matthysen | Multicausality of dispersal: a review, in Dispersal ecology and evolution[END_REF]. This intra-specific variation of behavioural responses may influence a wide range of eco-evolutionary processes [START_REF] Sih | Ecological implications of behavioural syndromes[END_REF][START_REF] Ronce | Dispersal syndromes, in Dispersal ecology and evolution[END_REF][START_REF] Edelaar | Non-random gene flow: an underappreciated force in evolution and ecology[END_REF]. Inter-individual differences in habitat choice behaviours may be genotype-dependent and related to differential performance in specific habitats or niches [START_REF] Sih | Ecological implications of behavioural syndromes[END_REF][START_REF] Edelaar | Non-random gene flow: an underappreciated force in evolution and ecology[END_REF][START_REF] Hawthorne | Genetic linkage of ecological specialization and reproductive isolation in pea aphids[END_REF][START_REF] Cousyn | Rapid, local adaptation of zooplankton behavior to changes in predation pressure in the absence of neutral genetic changes[END_REF][START_REF] De Meester | Evolutionary potential and local genetic differentiation in a phenotypically plastic trait of a cyclical parthenogen, Daphnia magna[END_REF][START_REF] Jaenike | Genetic variation for habitat preference: evidence and explanations[END_REF]. Genetic variation in habitat choice behaviour is common in animals and has been well documented for host plant preferences and cryptic habitat selection in insects (reviewed in [START_REF] Jaenike | Genetic variation for habitat preference: evidence and explanations[END_REF].

One of the cues animals use to select their habitat is light which animals may be positively or negatively attracted to. For example, negative phototaxis acts as a predator-avoidance mechanism in some aquatic organisms [START_REF] De Meester | Evolutionary potential and local genetic differentiation in a phenotypically plastic trait of a cyclical parthenogen, Daphnia magna[END_REF][START_REF] Cousyn | Rapid, local adaptation of zooplankton behavior to changes in predation pressure in the absence of neutral genetic changes[END_REF][START_REF] Michels | Inter-clonal variation in phototactic behaviour and key lifehistory traits in a metapopulation of the cyclical parthenogen Daphnia ambigua: the effect of fish kairomones[END_REF][START_REF] Borowsky | Responses to light in two eyeless cave dwelling amphipods (Niphargus ictus and Niphargus frasassianus)[END_REF] and may help some soil organisms that are very sensitive to relative humidity to prevent desiccation by looking for deeper and more humid soil layers [START_REF] Salmon | Linking species, traits and habitat characteristics of Collembola at European scale[END_REF][START_REF] Salmon | Species traits and habitats in springtail communities: A regional scale study[END_REF]. Nevertheless, different individuals or populations within the same species may present different phototactic preferences. For example, clonal populations of Daphnia magna (a zooplancton species) exposed to high levels of predatory pressure are more photophobic than clonal populations less exposed to predation [START_REF] Cousyn | Rapid, local adaptation of zooplankton behavior to changes in predation pressure in the absence of neutral genetic changes[END_REF][START_REF] De Meester | Evolutionary potential and local genetic differentiation in a phenotypically plastic trait of a cyclical parthenogen, Daphnia magna[END_REF]. In general, differences in phototaxis may have a heritable, genetic basis (e.g., [START_REF] Markow | Genetic analysis of phototactic behavior in Drosophila simulans[END_REF][START_REF] De Meester | Evolutionary potential and local genetic differentiation in a phenotypically plastic trait of a cyclical parthenogen, Daphnia magna[END_REF][START_REF] Cousyn | Rapid, local adaptation of zooplankton behavior to changes in predation pressure in the absence of neutral genetic changes[END_REF] or could be the result of non-genetic, phenotypic plasticity and personality (is this too anthropomorphic for Daphnia or springtails?) differences among individuals (e.g., [START_REF] Kain | Phototactic personality in fruit flies and its suppression by serotonin and white[END_REF]. Quantifying sources of variation in phototaxis is therefore important to understand the evolution of this widespread behavioural trait. Here, we investigate the genetic variation in light sensitivity and light-dependent movement behaviour in the eye-less springtail Folsomia candida Willem 1902 (Collembola, Isotomidae), an hemi-edaphic and cosmopolitan soil organism inhabiting various habitats such as caves, forest litter and man-made habitats [START_REF] Fountain | Folsomia candida (Collembola): a "standard" soil arthropod[END_REF].

The degeneration or even loss of the visual system is a convergent and frequent evolutionary phenomenon in soil-dwelling and cave animals [START_REF] Christiansen | Morphological adaptations, in Encyclopedia of caves[END_REF]. Nevertheless, even eyeless and eye-reduced species often retain some sensitivity to the ambient light level through extra-ocular photoreceptors (EOP, Taddei-Ferretti and Musio, 2000;[START_REF] Ullrich-Lüter | Unique system of photoreceptors in sea urchin tube feet[END_REF], which are useful for the maintenance of circadian rhythms [START_REF] Friedrich | Biological clocks and visual systems in cave-adapted animals at the dawn of speleogenomics[END_REF] or for orientation and habitat choice [START_REF] Timmermann | Phototactic response and light sensitivity in an epigean and a hypogean population of a barb (Garra barreimiae, Cyprinidae)[END_REF][START_REF] Borowsky | Responses to light in two eyeless cave dwelling amphipods (Niphargus ictus and Niphargus frasassianus)[END_REF]. Indeed, previous works strongly suggest that F. candida is sensitive to light despite being eyeless. In choice-test experiments, F. candida avoids UV light moving to warmer locations exposed to white light, prefers darkness over cool white light [START_REF] Fox | Phototactic responses to ultraviolet and white light in various species of Collembola, including the eyeless species, Folsomia candida[END_REF], and displays a dose-response avoidance of UV-B light relative to darkness [START_REF] Beresford | Lethal and sub-lethal effects of UV-B radiation exposure on the collembolan Folsomia candida (Willem) in the laboratory[END_REF]). Yet, to our knowledge, no study has examined the wavelengths of maximum sensitivity of the ocular or extra-ocular photoreceptors (EOPs)of these animals [START_REF] Barra | Les photorecepteurs des Collemboles, etude ultrastructurale[END_REF][START_REF] Jordana | A new type of arthropod photoreceptor[END_REF][START_REF] Fox | Phototactic responses to ultraviolet and white light in various species of Collembola, including the eyeless species, Folsomia candida[END_REF]. In true insects (Pterygota), few species are able to detect wavelengths longer than 600 nm (red light), which suggests a red-blind common ancestor [START_REF] Briscoe | The evolution of color vision in insects[END_REF]. In addition, the behavioural tests mentioned above could not always prevent confounding effects of differences in temperature or humidity associated with the lighting treatment. This is of great importance, as F. candida needs a relative humidity close to saturation [START_REF] Holmstrup | Strategies for cold and drought tolerance in permeable soil invertebrates, Ministry of Environment[END_REF][START_REF] Waagner | Recovery of reproduction after drought in the soil living Folsomia candida (Collembola)[END_REF] and is very sensitive to temperature (Boiteau andMackinley, 2012, 2013).

Although sexual reproduction exists in some populations of F. candida [START_REF] Frati | High levels of genetic differentiation between Wolbachia-infected and non-infected populations of Folsomia candida (Collembola, Isotomidae)[END_REF], this species is generally recognized as asexual, and most studies using F. candida as a model species have used parthenogenetic lineages [START_REF] Fountain | Folsomia candida (Collembola): a "standard" soil arthropod[END_REF]. Earlier studies on several parthenogenetic lineages have uncovered substantial intra-specific genetic and morphological polymorphism [START_REF] Chenon | Genetic polymorphism in nine clones of a parthenogenetic collembolan used in ecotoxicological testing[END_REF][START_REF] Tully | Two major evolutionary lineages revealed by molecular phylogeny in the parthenogenetic collembola species Folsomia candida[END_REF][START_REF] Tully | Intraspecific phenotypic variation and morphological divergence of strains of Folsomia candida (Willem) (Collembola: Isotomidae), the "standard" test springtail[END_REF]. Intra-specific diversity is organised in two major evolutionary clades [START_REF] Tully | Two major evolutionary lineages revealed by molecular phylogeny in the parthenogenetic collembola species Folsomia candida[END_REF][START_REF] Tully | Intraspecific phenotypic variation and morphological divergence of strains of Folsomia candida (Willem) (Collembola: Isotomidae), the "standard" test springtail[END_REF], and life history studies have shown that two contrasted biodemographic strategies evolved along the divergence of these two clades [START_REF] Tully | Reproductive flexibility: genetic variation, genetic costs and long-term evolution in a collembolan[END_REF][START_REF] Tully | Facteurs génétiques, maternels et environnementaux de l'expression des traits d'histoire de vie chez le collembole[END_REF][START_REF] Tully | The evolution of postreproductive life span as an insurance against indeterminacy[END_REF][START_REF] Mallard | Within species variation in long-term trajectories of growth, fecundity and mortality in the Collembola Folsomia candida[END_REF]. One clade has a high reproductive potential: when sufficient food is available, these springtails produce on average larger clutches than the ones from the other clade [START_REF] Tully | Reproductive flexibility: genetic variation, genetic costs and long-term evolution in a collembolan[END_REF], but they have shorter mean lifespans than the less fecund clade and also reach a smaller adult body size [START_REF] Tully | Reproductive flexibility: genetic variation, genetic costs and long-term evolution in a collembolan[END_REF][START_REF] Tully | The evolution of postreproductive life span as an insurance against indeterminacy[END_REF][START_REF] Mallard | Within species variation in long-term trajectories of growth, fecundity and mortality in the Collembola Folsomia candida[END_REF]. These two groups of clonal lineages fit well to the typical slow (A) and fast (B) life history syndromes (or r-K life histories, see [START_REF] Reznick | R-and K-selection revisited: the role of population regulation in life-history evolution[END_REF]. But, until now, the ecological conditions in which they have evolved and the time elapsed since the divergence of the two clades remain to be determined. Intra-specific variation in habitat choice behaviour and mobility has so far neither been examined in this species nor in other Collembola. Instead, the few works that relate the habitat preference or distribution and colonization ability of Collembola with their morphological and life history traits are focused on the study of collembolan community composition [START_REF] Salmon | Linking species, traits and habitat characteristics of Collembola at European scale[END_REF]Ponge andSalmon, 2013, Huebner et al, 2012;[START_REF] Salmon | Species traits and habitats in springtail communities: A regional scale study[END_REF][START_REF] Ponge | Decreased biodiversity in soil springtail communities: the importance of dispersal and landuse history in heterogeneous landscapes[END_REF].

Intraspecific variation in phototactic behaviour and life history traits has been well investigated in Daphnia magna. In this species, positive phototactic clones present a fast life history strategy whereas negative and intermediate phototactic clones present a slow life history strategy (e.g., De Meester, 1994).

We tested if springtail clonal variation in light-dependent habitat choice exists using an experimental setup to control the lighting conditions while maintaining constant temperature and moisture. We tracked springtail movements under this setup to quantify their spatial preference for shaded versus illuminated areas as a measure of the habitat choice behaviour. We first tested whether F. candida can use light as an environmental cue for habitat choice under different lighting conditions, including natural shaded and sunny understory spectra and an artificial red-light spectrum. We measured springtail sensitivity to long wavelengths because these wavelengths are dominant under the forest canopy [START_REF] Smith | Sensing the light environment: the functions of the phytochrome family[END_REF] while red and far red are the principal wavelengths that penetrate the soil [START_REF] Bliss | Penetration of light into soil and its role in the control of seed germination[END_REF]. We hypothesised that F. candida should not detect or react to red and far-red light per se if the incapability to detect these wavelengths was a shared condition of most true insects and Collembola, even though most springtails are sensitive to heat generated by red light [START_REF] Briscoe | The evolution of color vision in insects[END_REF]. We further studied the sensitivity to light of twelve clonal lineages of F. candida including eleven lineages belonging to the two evolutionary clades described earlier [START_REF] Tully | Two major evolutionary lineages revealed by molecular phylogeny in the parthenogenetic collembola species Folsomia candida[END_REF][START_REF] Tully | Reproductive flexibility: genetic variation, genetic costs and long-term evolution in a collembolan[END_REF][START_REF] Tully | The evolution of postreproductive life span as an insurance against indeterminacy[END_REF][START_REF] Tully | Intraspecific phenotypic variation and morphological divergence of strains of Folsomia candida (Willem) (Collembola: Isotomidae), the "standard" test springtail[END_REF]. We addressed the following questions: Do light sensitivity and habitat choice behaviour vary between clonal lineages, as has been found in other taxa [START_REF] Jaenike | Genetic variation for habitat preference: evidence and explanations[END_REF][START_REF] De Meester | Evolutionary potential and local genetic differentiation in a phenotypically plastic trait of a cyclical parthenogen, Daphnia magna[END_REF][START_REF] Cousyn | Rapid, local adaptation of zooplankton behavior to changes in predation pressure in the absence of neutral genetic changes[END_REF])? If such clonal variation exists, how is it organized relative to the phylogenetic clades and what are the links between the behavioural responses and the main life history strategies of each clade? We predicted that lineages from the slow life history group would be more photophobic [START_REF] De Meester | Life histories and habitat selection in Daphnia: divergent life histories of D. magna clones differing in phototactic behavior[END_REF][START_REF] De Meester | Life history characteristics of Daphnia magna clones differing in phototactic behavior[END_REF], given that photophobia is likely to be associated with life in more stable habitats, which usually selects for a slow life history [START_REF] Pianka | R-Selection and K-Selection[END_REF][START_REF] Reznick | R-and K-selection revisited: the role of population regulation in life-history evolution[END_REF].

Materials and methods

Maintenance and origin of the studied springtails

We used twelve clonal lineages of the Collembola Folsomia candida labelled AP, BR, BV, DK, GB, GM, HA, ME, PB, TO, US and WI [START_REF] Tully | Two major evolutionary lineages revealed by molecular phylogeny in the parthenogenetic collembola species Folsomia candida[END_REF]. Information about the phylogenetic relationships and habitat and geographical origin of all strains except ME can be found in previous studies [START_REF] Tully | Two major evolutionary lineages revealed by molecular phylogeny in the parthenogenetic collembola species Folsomia candida[END_REF][START_REF] Tully | Intraspecific phenotypic variation and morphological divergence of strains of Folsomia candida (Willem) (Collembola: Isotomidae), the "standard" test springtail[END_REF]. Clones AP, BV, BR, GB and HA belong to the "slow clade" A while the "fast clade" B comprises the clones DK, GM, PB, TO, US and WI [START_REF] Tully | Reproductive flexibility: genetic variation, genetic costs and long-term evolution in a collembolan[END_REF]. The new clone ME was collected in November 2013 from some decaying wood beams into an abandoned man-made tunnel in the Mercantour French National Park (South-East of France, 44° 7.026'N, 7° 16.727'E, 1530 m). The life history strategy of this clone and its phylogenetic relationships with the other clones are currently unknown.

All clonal populations were reared in similar conditions in polyethylene vials (inner diameter 52 mm, height 65 mm) filled with a 30 mm layer of plaster of Paris mixed with Indian ink to increase visual detectability of individuals [START_REF] Tully | Reproductive flexibility: genetic variation, genetic costs and long-term evolution in a collembolan[END_REF]. Populations were kept in incubators at 20°C (± 0.5°C) in the dark and fed with pellets of a mixture of agar and dried yeast [START_REF] Tully | Reproductive flexibility: genetic variation, genetic costs and long-term evolution in a collembolan[END_REF]. We established synchronised populations of each clone in the same way and at the same time by transferring 10-12 randomly chosen adult females from stock cultures to new culture vials. Females were transferred to new vials every week and old vials were kept at 20°C for laid eggs to hatch since this temperature is optimal for all lineages. Then, vials were checked weekly for the presence of new-borns. If new-borns were detected, ad libitum food was provided weekly to ensure optimum body growth. Therefore, the age of the experimental animals was known with a one-week accuracy. However, due to the low fecundity of some clones or accidental flooding in some vials, we had to add medium size (~1.4mm) adults of unknown age taken from the stock cultures in order to keep balanced samples for all clones and treatments (see Table S1 for details on sample size and individual characteristics).

Measurement of natural light conditions

At midday of a sunny day in June 2013, we measured 23 light spectra under the canopy of a temperate forest located in the Foljuif field station, 80km South of Paris, France (48°17'13.96"N, 2°40'40.34"E), where populations of F. candida have been previously found. We recorded 19 spectra in the shade and four spectra on sunspots in the undergrowth. The climate conditions and plant community structure of this forest have been studied extensively [START_REF] Blandin | Etude d'un ecosysteme forestier mixte, 1: Climat, structure de la vegetation et retombees de litieres [France, Seine-et-Marne[END_REF]. We used a handheld spectrometer parameterized for absolute irradiance measurements in the range 200-850 nm (Jaz Series, JAZ-ULM-200, Ocean Optics, USA). We measured the incident light with a cosine-corrected probe (180° field of view) after dark calibration and converted data into units of irradiance spectrum (µW.cm -2 ) using factory calibration of the spectrophotometer. Measurements were taken at ground levels under the canopy of different tree and bush species, as well as under sunspots under the trees and into tree stump holes to sample different light spectra that natural populations of springtails are likely to encounter aboveground.

Simulation of environmental conditions

Artificial light and constant climate conditions were simulated in a controlled environment laboratory at the CEREEP EcotronIleDeFrance (Saint-Pierre-les-Nemours, France). We used one 13 m 3 controlled environment chamber of the Ecolab, where climate (temperature, relative humidity, rainfall) and lighting conditions can be simulated [START_REF] Verdier | Climate and Atmosphere Simulator for Experiments on Ecological Systems in Changing Environments[END_REF]. In the centre of the chamber, we installed vials for the characterisation of photophobia in springtails (see below) on a controlled water table to maintain constant humidity and temperature conditions in the vials during the observations. Controlling for constant humidity and temperature was essential to test exclusively for the effect of light, since slight changes in humidity and temperature could modify the springtails' behaviour and thus affect their apparent or real photophobia [START_REF] Holmstrup | Strategies for cold and drought tolerance in permeable soil invertebrates, Ministry of Environment[END_REF][START_REF] Boiteau | Locomotor response of Folsomia Candida (Collembola: Isotomidae) to cooling temperatures[END_REF]. Atmospheric temperature in the chamber was maintained at 20ºC (± 0.1°C) and temperature of the water table was set at 10ºC. We used a cold water table to keep the temperature inside the vials slightly lower than the air temperature to prevent condensation on the transparent lid that covered the vials (Figure 1). The temperature measured inside the vials was 19ºC and the relative humidity was 93% (DS-1923 iButton loggers) irrespective of the test conditions.

Light was provided by means of modular LED (light emitting diodes) arrays (n=40) allowing to turn on or off cool white LEDs (n=15 per array) and five different LED types with maximum emission wavelengths in the UV (370 nm, n=3), green (520 nm, n=2), red (660 nm, n=5), far red (740 nm, n=7), and infrared (840 nm, n=3). Thus, up to 36,864 combinations of LEDs parameters can be programmed. To automatize the simulation of the most appropriate lighting conditions, we wrote a procedure in R (R Development Core Team 2012) to calculate the sum of squared distance between the reference spectrum and the simulated spectrum, and to select the combination of parameters that minimised this sum. We simulated the "average undergrowth" light spectrum by producing a spectrum which best matched the mean spectra of the 23 light measurements made at various places on the forest floor (see Figure 2a and Table S2 provided as supplementary information). We also produced a "maximum undergrowth spectrum", using the four measurements made on sunspots, which should correspond to the strongest undergrowth light conditions that litterdwelling springtails are susceptible to be exposed to in this forest (Figure 2b and Table S2). We expected a stronger light-avoidance response under this condition compared to under the mean undergrowth light spectrum. In addition, we simulated a "red light" spectrum whose power was close to that of the mean undergrowth light spectrum (Figure 2a and Table S2). This red light spectrum has two irradiance peaks (red at 660 nm and far red at 740 nm) and has no UV light that could promote photophobia in springtails [START_REF] Fox | Phototactic responses to ultraviolet and white light in various species of Collembola, including the eyeless species, Folsomia candida[END_REF][START_REF] Beresford | Lethal and sub-lethal effects of UV-B radiation exposure on the collembolan Folsomia candida (Willem) in the laboratory[END_REF]. We further compared these three lighting conditions with a control treatment where the room was maintained in the dark. This condition was used to verify whether the experimental set-up in itself could affect the behaviour of the springtails and to compare the mobility of the different clones in the dark.

Experimental protocol

The dark, average undergrowth and red light spectra were tested in June 2014, and the maximum undergrowth spectrum was tested in December 2014. We used 36 individuals of each clonal lineage in each light treatment (n = 1728), and six individuals per clone were tested every day. Springtails were photographed to measure their length (mean = 1.41 mm, range =0.93-2.19 mm) and kept isolated without food in rearing vials the day before. These individuals were transferred the next morning to new rearing vials wrapped with a thin flexible black plastic sheet, and only exposed to light from above. Each plastic box was filled with plaster of Paris and had a hole in the bottom to ensure that the plaster can be kept constantly wet to saturation (Figure 1). An opaque black piece of plastic covered half of the vial's opening and descended vertically along the centre of the vial down to 3 mm above the plaster level, in order to maintain in the dark half of the vial while allowing the springtail to move freely between the two halves of the vials. Preliminary experiments showed that in total darkness springtail behaviour is very sensitive to drought and slight humidity differences between the two sides of the vials. The vials were thus covered with transparent lids to prevent any desiccation. When being transferred, animals were separated randomly in two groups of same size and released in the bright and dark sides at the beginning of the experiment (here called 'position') respectively, such that the proportion of animals in each position was exactly 50% for each clone and treatment. One observer (MGR) studied the springtails' movements between the two sides of the vials by recording the position of each animal every hour during six hours (from 9 am to 5 pm).

For the control dark treatment, we observed the springtails' positions as quickly as possible using a soft light head torch to minimise disturbance. At the end of the experiment, we inspected all vials to check that no individual had escaped or deceased to ensure valid observations when the springtails did not move.

Statistical analyses

We used the software R (R Development Core Team 2012) for the statistical analysis and to produce the graphs. The raw data and the R-scripts of the analysis, and figures are provided as supplementary information. We calculated the movement (or transition) probabilities (the probability that an individual has moved from one side to the other side of the vial after one hour), using pairs of consecutive observations of the same individual to form a binomial variable which equals zero when the individual stayed in the same half of the vial and equals one when the animal had moved to the other. This binomial variable was analysed with nested generalised linear mixed models (GLMMs) for repeated measurements approach with a binomial error distribution and logit link function, using treatment (light spectrum type), position (dark or light) and clonal group as the fixed factors, and individual as the random factor. Observations were thus clustered in individuals and individuals were nested within clones. In the full model, we analysed all main effects as well as two and three-way interactions. We implemented the full model with the glmer procedure in package lme4 for R 3.1.2 (Bates et al, 2012; R Development Core Team 2012) and tested for the significance of fixed effects with Chi-square tests. Since the three-way interaction Spectrum*Position*Clone was significant (see Table S3 provided as supplementary information), we analysed each treatment independently to get a better understanding of the behavioural differences among clones under each treatment.

For subsequent analyses we further included the body length of each individual in the analysis to test for a size effect on individual movement probability. Since mean body length can vary substantially among clones [START_REF] Mallard | Within species variation in long-term trajectories of growth, fecundity and mortality in the Collembola Folsomia candida[END_REF], we calculated for each treatment the corrected body lengths as the sum of the overall mean body length and the body length residuals computed here as the difference between individual lengths and mean length of their clone. This "corrected body length" was then used to study the effect of body length while controlling for differences among clones. The full models included position, clone and corrected length as fixed factors, the two-way interactions and individual identity as the random factor. We used a backward stepwise selection to get final models that only included significant factors and interactions. We also report full models in Table S4 following Forstmeier and Schielzeth's recommendation [START_REF] Forstmeier | Cryptic multiple hypotheses testing in linear models: overestimated effect sizes and the winner's curse[END_REF]. For the full and final models, we calculated marginal and conditional R 2 , which quantify the amount of variance explained by the models [START_REF] Nakagawa | A general and simple method for obtaining R2 from generalized linear mixed-effects models[END_REF]. We additionally performed GLMMs with clade, position and their interaction as fixed factors, and clone and individual within clone as nested random factors to test for phylogenetic variation. The clone ME, whose evolutionary clade has not yet been characterized, was excluded from these analyses.

Results

In the homogenous dark environment, the mobility of springtails was not influenced by their position in the vial (Table 1; Figure 3) and was independent of their body length. The individual mobility (average transition probability) differed significantly among clones (Table 1, Figure 4a), and this variability was mostly due to difference between the two evolutionary clades (χ 2 1= 15.86, p<0.001). In complete darkness, the clones with a fast life history (clade B) had on average a higher propensity to switch between sides (0.44 per hour) than clones with a slow life-history strategy (0.34, Fig 4a).

Contrary to our expectations, F. candida proved to be slightly sensitive to red light (Table 2), independently of their length: springtails preferentially moved from red light to dark (0.46 per hour) than from dark to red light (0.42, see Figure 3). Clones also differed in their mean mobility, but these differences were not due to clade (χ 2 1 =0.07, p = 0.79, Figure 4b).

Springtails showed photophobia when exposed to average undergrowth light: they moved from the light to the dark side significantly more often than in the reverse direction (0.43 versus 0.37, see Fig 3 andTable 3). The different clones (χ 2 11= 18.09, p=0.08) and the two clades (χ 2 1= 0.40, p=0.40) had similar movement probabilities (Figure 4c and Table 3). In this light condition, movement probability varied with body length depending on position: while in the dark compartment the length did not influence mobility, but when exposed to light, smaller than average springtails were more prone to move than longer ones (Figure 5a). When exposed to maximum undergrowth light, springtails moved away on average more often (0.45 per hour) than when they were in the dark side of the vials (0.36, see Figure 3 and Table 4). Differences among clones (χ 2 11=15.67, p=0.15) and clades were not significant (χ 2 1 =1.04, p=0.31, Figure 4d), and there was a non-significant trend (Position* Corrected Length: χ 2 11=1.22, p=0.27) for longer than average individuals to move less often than smaller individuals, as observed in the average undergrowth light treatment (Table 4). However, there was a significant interaction between effects of clone and of corrected length: for most clones, size matters little and smaller individuals tended to move more, while for clone GB the mobility was positively influenced by relative length (χ 2 11=21.49, p=0.03) (Figure 5b)

Clone GB was characterized by the slowest life history strategy.

Discussion

Limits of the experimental protocol

Our experimental protocol allowed us to study the behaviour of isolated individuals in darkness and under three light conditions while maintaining constant temperature and humidity. As the springtails' mobility is very sensitive to these environmental factors, this represents a significant improvement to test for photophobia.

Differences in temperature between the illuminated and dark sides of the boxes fell within the measurement error of our thermometer devices (< 0.5°C). Thus, even though we cannot strictly exclude that the bright side may be slightly warmer than the dark side, the potential temperature difference was so small that we could not detect it. Thus we interpret our result as a light sensitivity but one has to remember that temperature gradients could also play a role. For instance, temperature receptors such as micro setae on the antennas or legs may help springtails to perceive small temperature gradients, in particular when doing 'turning' movements (loops), a frequent behaviour when in an unknown environment [START_REF] Bengtsson | Irregular walks and loops combines in small-scale movement of a soil insect: implications for dispersal biology[END_REF],

In addition, our experimental design is not suited to know how many times springtails moved between sides during the one hour census interval. Our estimations of transition probabilities are based on the hypothesis that each springtail does no or only one move in an hour, which may be wrong given that some individuals may have moved multiple times between two censuses. Future studies should track more continuously individuals through space to obtain better estimates of dispersal distances and capacities.

Light sensitivity and its ecological significance

We found that, despite being eyeless, F. candida is slightly but significantly sensitive to light.

Photophobia was detected in the presence of the three light sources, including two light spectra chosen to mimic natural conditions encountered by springtails. Hence, this study confirms earlier general findings of negative phototaxis in this species [START_REF] Fox | Phototactic responses to ultraviolet and white light in various species of Collembola, including the eyeless species, Folsomia candida[END_REF][START_REF] Beresford | Lethal and sub-lethal effects of UV-B radiation exposure on the collembolan Folsomia candida (Willem) in the laboratory[END_REF][START_REF] Boiteau | Geotaxis in the springtail Folsomia candida[END_REF][START_REF] Salmon | Responses to light in a soil-dwelling springtail[END_REF]. But our study goes one step further by

showing that the intensity of photophobia response varied depending on the light spectrum: Light avoidance by springtails increased from the red light spectrum at an intensity equivalent to that of the average light spectrum under the canopy, to average undergrowth spectrum, and to maximum undergrowth spectrum. F. candida is thus sensitive both to light intensity and to its quality. This suggests that despite being blind, these springtails can use light as an environmental cue for orientation, movement and habitat choice under natural conditions. Light avoidance behaviour was not very strong, because individuals also moved from the dark to the illuminated side of the vials and switched their position frequently. This may be due to the fact that the intensity of the light spectrum tested was low compared for example to direct sunlight, but could also be related to the small spatial scale of the setup and the short time scale of our tests.

Multiple examples of eye-reduced and eyeless species that retain brightness and colour sensitivity exist (e.g. [START_REF] Ramirez | Understanding the dermal light sense in the context of integrative photoreceptor cell biology[END_REF][START_REF] Friedrich | Biological clocks and visual systems in cave-adapted animals at the dawn of speleogenomics[END_REF]. Extra-ocular photoreceptors (EOP) or non-visual photoreceptors can be responsible of light sensitivity when eyes are absent. EOP exist in the form of dermal photoreceptors in invertebrates [START_REF] Tosini | Dermal photoreceptors regulate basking behavior in the lizard Podarcis muralis[END_REF][START_REF] Binder | The role of dermal photoreceptors during the sea lamprey (Petromyzon marinus) spawning migration[END_REF][START_REF] Xiang | Light-avoidance-mediating photoreceptors tile the Drosophila larval body wall[END_REF][START_REF] Ullrich-Lüter | Unique system of photoreceptors in sea urchin tube feet[END_REF], the sixth caudal ganglion in decapods [START_REF] Wilkens | Photosensitivity in the sixth abdominal ganglion of decapod crustaceans: a comparative study[END_REF][START_REF] Larimer | A functional caudal photoreceptor in blind cavernicolous crayfish[END_REF], genitalia [START_REF] Arikawa | Multiple extraocular photoreceptive areas on genitalia of butterfly Papilio xuthus[END_REF] and neurons of the optic lobes in butterflies [START_REF] Lampel | Expression of UV-, blue-, long-wavelengthsensitive opsins and melatonin in extraretinal photoreceptors of the optic lobes of hawk moths[END_REF], and photosensitive cells in the hydra [START_REF] Taddei-Ferretti | Photobehaviour of Hydra (Cnidaria, Hydrozoa) and correlated mechanisms: a case of extraocular photosensitivity[END_REF]. F. candida lacks external eye facets and any other external eye structures as revealed by scanning electron microscopy [START_REF] Fox | Phototactic responses to ultraviolet and white light in various species of Collembola, including the eyeless species, Folsomia candida[END_REF]. Whether F. candida is only facet-less but has internal photoreceptors where eyes should be or whether it uses other types of extra-ocular photoreception such as a dispersed photoreception dermal sense [START_REF] Ramirez | Understanding the dermal light sense in the context of integrative photoreceptor cell biology[END_REF] or neural photoreception [START_REF] Lampel | Expression of UV-, blue-, long-wavelengthsensitive opsins and melatonin in extraretinal photoreceptors of the optic lobes of hawk moths[END_REF][START_REF] Wilkens | Photosensitivity in the sixth abdominal ganglion of decapod crustaceans: a comparative study[END_REF], remains an open question.

Photophobia may help springtails to avoid open spaces when they forage in the litter close to the surface, and together with positive geotaxis [START_REF] Boiteau | Geotaxis in the springtail Folsomia candida[END_REF], it may serve for orientation in vertical movements. Being able to avoid open spaces and the surface of the litter may provide several benefits. First of all, all clonal lineages tested here were non-pigmented and thus little protected against the deleterious effects of UV radiations. Long-term (two weeks) exposure to UV-B light increases mortality and causes DNA damage in this species [START_REF] Beresford | Lethal and sub-lethal effects of UV-B radiation exposure on the collembolan Folsomia candida (Willem) in the laboratory[END_REF].

Behavioural tactics of habitat selection for dark sheltered zones or deeper soil layers can help springtails to escape from potentially dangerous UV-radiation. This may also explain why unpigmented or slightly pigmented Collembola are preferentially distributed in edaphic or hemiedaphic rather than epigeal habitats and in woodlands rather than in grasslands [START_REF] Salmon | Linking species, traits and habitat characteristics of Collembola at European scale[END_REF][START_REF] Salmon | Species traits and habitats in springtail communities: A regional scale study[END_REF]. However, the photophobia of F. candida cannot be exclusively related to UV-radiation avoidance, since F. candida was also sensitive to red light. Long-waves are dominant wavelengths under the forest canopy [START_REF] Smith | Sensing the light environment: the functions of the phytochrome family[END_REF], and thus sensitivity to these wavelengths could be related to a more general aboveground avoidance. Moreover, as ambient light level is often correlated with other climatic factors (notably temperature and humidity), photophobia is useful to target deeper, more humid and more stable microhabitats preferred by most Collembola species, especially those that are the most vulnerable to desiccation [START_REF] Salmon | Linking species, traits and habitat characteristics of Collembola at European scale[END_REF][START_REF] Salmon | Species traits and habitats in springtail communities: A regional scale study[END_REF]). Another explanation is that photophobia may indirectly help springtails to avoid predators.

Being devoid of compound eyes F. candida probably cannot form images and thus cannot see approaching predators [START_REF] Hauzy | Spatial heterogeneity and functional response: an experiment in microcosms with varying obstacle densities[END_REF], as it is the case of the sister eyeless species Folsomia fimentaria [START_REF] Baatrup | Predation of the mite Hypoaspis aculeifer on the springtail Folsomia fimetaria and the influence of sex, size, starvation, and poisoning[END_REF]. This species is thus very vulnerable in illuminated environments to visually active-chasing predators such as hunting ground-beetles, some of which are ferocious predators of springtails [START_REF] Ernsting | Effects of Food Deprivation and Type of Prey on Predation by Notiophilus biguttatus F. (Carabidae) on Springtails (Collembola)[END_REF].

Intra-specific variability in movement and light sensitivity

In general, smaller individuals were more mobile than larger ones and we found that smaller springtails exposed to light were also more likely to seek refuge in the dark side than the larger ones in the average undergrowth light treatment. The risk of staying in an illuminated environment may not be the same for large and small individuals. For example, it has been shown that the carabid Notiophilus biguttatusa common visual hunter and voracious predator of springtails -preys more efficiently on small rather large Collembola [START_REF] Ernsting | Components of Predatory Behaviour Underlying Density-Dependent Prey-Size Selection by Notiophilus biguttatus F. (Carabidae, Coleoptera)[END_REF]. This interpretation is tenuous given that the effect of body size on movement was small, variable across clones and not repeatable across light treatments. These variable effects of body size are difficult to explain. The intensity of photophobia may also depend on other individual state parameters such as the nutritional status or the phase in the moulting and oviposition cycle. There could also be a trade-off between the preference for safer dark deep zones and the need for foraging and dispersal at the illuminated surface [START_REF] Dromph | Effect of starvation on phototaxis and geotaxis of collembolans[END_REF].

We searched for clonal variation in both the mean movement propensity and in the light sensitivity of movement. We first found that the various lineages moved differently in the dark: the average transition probabilities varied from 0.2 to almost 0.5 depending on the lineages. Quite remarkably, this genetic diversity was structured between the two previously described phylogenetic clades: clones from the clade characterised by a fast life-history strategy had on average higher transition probabilities than those with slow life-history strategy, suggesting that they move faster and can potentially disperse over longer distances in natural conditions. Even though we did not measure dispersal per se (defined as the movement of individuals or propagules that leads to gene flow, see [START_REF] Clobert | Dispersal ecology and evolution, 1st ed[END_REF], mobility or exploratory activity has generally been positively correlated to dispersiveness and can thus be used as a rough proxy for dispersal in natural populations [START_REF] Ronce | Dispersal syndromes, in Dispersal ecology and evolution[END_REF]. The higher movement propensity of the fast life-history clade agrees with the evolutionary scenario of the 'colonizer syndrome' where fast life-history species or populations are more dispersed than their slow life-history counterparts [START_REF] Ronce | Dispersal syndromes, in Dispersal ecology and evolution[END_REF]. The combination of increased dispersal, high fecundity but short lifespan and low competitive ability allows fast individuals to colonise new habitats more efficiently. Individuals with reduced movement propensity, low fecundity high competitive ability perform better at the core range and in stable habitats, where populations have reached their carrying capacity and where intraspecific competition is strong [START_REF] Burton | Trade-offs and the evolution of life-histories during range expansion[END_REF]. Future studies should assess whether differences in short-distance movement probabilities translate into clonal differences in longdistance movement behaviour and therefore dispersal distances.

Studies on the intraspecific variation of phototaxis are scarce but differences in phototaxis among genetically distinct lineages have been reported in several species. It has been shown that several Daphnia clonal populations differ in their photophobia in response to contrasting predatory pressures [START_REF] Cousyn | Rapid, local adaptation of zooplankton behavior to changes in predation pressure in the absence of neutral genetic changes[END_REF][START_REF] De Meester | Evolutionary potential and local genetic differentiation in a phenotypically plastic trait of a cyclical parthenogen, Daphnia magna[END_REF]. Artificially selected Drosphila strains also differ in their phototactic preferences both in quality and degree as a result of different genotypes and the presence of specific mutations [START_REF] Kain | Phototactic personality in fruit flies and its suppression by serotonin and white[END_REF]. But contrary to these findings, we did not find any significant differences in our model organism in light sensitivity among clones and clades for the two undergrowth spectra, including when springtails were exposed to the brightest light spectrum.

This could indicate that photophobia is a shared, evolutionary rigid behavioural trait in this species.

It also implies that other studies that failed to find intraspecific variation in phototaxis may be missing in the published bibliography.

In addition, when exposed to light, the genetic differences in movement behaviour observed in the dark vanish. Our observations also suggest that moving from one side of the vial to the other side in homogeneous dark conditions or in heterogeneous light conditions is associated with different underlying behaviours: exploration for the former and escape for the latter. This is also supported by the fact that body length did not affect movement in the dark while small individuals were found to be slightly more mobile than larger ones in the heterogeneous environments.

Figure 1 .

 1 Figure 1. a) Photographs of the experimental boxes. Boxes were wrapped with a thin flexible black

Fig 2 .Fig 3 .

 23 Fig 2. Light conditions used during the experiment. a) Natural and simulated average under-canopy

Fig 4 .

 4 Fig 4. Mean movement probability between the two sides of the test boxes (±95% CI) for each

Fig 5 .

 5 Fig 5. a) Predicted mean movement probability between the two sides as a function of individuals'

.

  Models that best describe transition probabilities in the average undergrowth spectrum after backward selection (N=432). Marginal and conditional R 2 of the final model are 0.011 and 0.074 respectively.

		Variance	Std. Dev. Chi squared	D.f. p-value
	Fixed factors:
	Position		9.6243	1	0.0019
	Residual length		3.8848	1	0.0487
	Position* Corrected length		6.9324	1	0.008
	Random factors:
	Individual	0.2255	0.4748
	Dropped variables:
	Clone* Corrected Length		12.631	11	0.3181
	Clone* Position		17.509	11	0.0937
	Clone		18.088	11	0.0796

Table 4 .

 4 Models that best describe transition probabilities in the maximum undergrowth spectrum after backward selection (N=432). Marginal and conditional R 2 of the final model are 0.039 and 0.089 respectively.

	Variance

Std. Dev. Chi squared D.f. p-value Fixed factors:

	Position		23.4988	1	<0.001
	Clone		15.6741	11	0.1536
	Corrected length		0.0312	1	0.8597
	Clone* Corrected length		21.4906	11	0.0286
	Random factors:
	Individual	0.1818	0.4263
	Dropped variables:
	Clone* Position		11.871	11	0.3734
	Position* Corrected Length		1.2208	11	0.2692

Acknowledgments

We want to thank to the staff members of the especially Simon Chollet,

Mathieu

Llavata and Florent Massol. This research was supported by the Centre National de la Recherche Scientifique and a grant from the Ecole doctorale Diversité du Vivant at the Université Pierre and Marie Curie to M.G.R. (contract 859/2013). This work was made possible thanks to the CNRS Research Infrastructure Ecotrons supported by the Conseil régional d'Ile-de-France (DIM R2DS I-05-098/R and 2011-11017735) and the Agence Nationale de la Recherche grant of the "Investissements d'avenir" program (ANR-11-INBS-0001 AnaEE France).


S Buitink

A Corstanje

H Falcke

J R Hörandel

T Huege

A Nelles

J P Rachen

L Rossetto

P Schellart

O Scholten

S Ter Veen

S Thoudam

T N G Trinh

J Anderson

A Asgekar

I M Avruch

M E Bell

M J Bentum

G Bernardi

P Best

A Bonafede

F Breitling

J W Broderick

W N Brouw

M Brüggen

H R Butcher

D Carbone

B Ciardi

J E Conway

F De Gasperin

E De Geus

A Deller

R.- J Dettmar

G Van Diepen

S Duscha

J Eislöffel

D Engels

J E Enriquez

R A Fallows

R Fender

C Ferrari

W Frieswijk

M A Garrett

J M Grießmeier

A W Gunst

M P Van Haarlem

T E Hassall

G Heald

J W T Hessels

M Hoeft

A Horneffer

M Iacobelli

H Intema

E Juette

A Karastergiou

V I Kondratiev

M Kramer

M Kuniyoshi

G Kuper

J Van Leeuwen

G M Loose

P Maat

G Mann

S Markoff

R Mcfadden

D Mckay-Bukowski

J P Mckean

M Mevius

D D Mulcahy

H Munk

M J Norden

E Orru

H Paas

M Pandey-Pommier

V N Pandey

M Pietka

R Pizzo

A G Polatidis

W Reich

H J A Röttgering

A M M Scaife

D J Schwarz

M Serylak

J Sluman

O Smirnov

B W Stappers

M Steinmetz

A Stewart

J Swinbank

M Tagger

Y Tang

C Tasse

M C Toribio

R Vermeulen

C Vocks

C Vogt

R J Van Weeren

R A M J Wijers

S J Wijnholds

M W Wise

O Wucknitz

S Yatawatta

P Zarka

J A Zensus

A large light-mass component of cosmic rays at 10 17 -10 17.5 eV from radio observations

Cosmic rays are the highest energy particles found in nature. Measurements of the mass composition of cosmic rays between 10 17 eV and 10 18 eV are essential to understand whether this energy range is dominated by Galactic or extragalactic sources. It has also been proposed that the astrophysical neutrino signal [START_REF] Aartsen | Evidence for High-Energy Extraterrestrial Neutrinos at the IceCube Detector[END_REF] comes from accelerators capable of producing cosmic rays of these energies [START_REF] Murase | Testing the hadronuclear origin of PeV neutrinos observed with IceCube[END_REF] . Cosmic rays initiate cascades of secondary particles (air showers) in the atmosphere and their masses are inferred from measurements of the atmospheric depth of the shower maximum, X max [START_REF] Aab | Depth of maximum of air-shower profiles at the Pierre Auger Observatory. II. Composition implications[END_REF] , or the composition of shower particles reaching the ground [START_REF] Apel | Ankle-like feature in the energy spectrum of light elements of cosmic rays observed with KASCADE-Grande[END_REF] .

Current measurements [START_REF] Kampert | Measurements of the cosmic ray composition with air shower experiments[END_REF] suffer from either low precision, and/or a low duty cycle. Radio detection of cosmic rays [START_REF] Allan | Radio Emission from Extensive Air Showers[END_REF][START_REF] Falcke | Detecting radio emission from cosmic ray air showers and neutrinos with a digital radio telescope[END_REF][START_REF] Falcke | Detection and imaging of atmospheric radio flashes from cosmic ray air showers[END_REF] is a rapidly developing technique [START_REF] Huege | The Renaissance of Radio Detection of Cosmic Rays[END_REF] , suitable for determination of X max 10, 11 with a duty cycle of in principle nearly 100%. The radiation is generated by the separation of relativistic charged particles in the geomagnetic field and a negative charge excess in the shower front [START_REF] Allan | Radio Emission from Extensive Air Showers[END_REF][START_REF] Werner | Macroscopic Treatment of Radio Emission from Cosmic Ray Air Showers based on Shower Simulations[END_REF] . Here we report radio measurements of X max with a mean precision of 16 g/cm 2 between 10 17 -10 17.5 eV. Because of the high resolution in X max we can determine the mass spectrum and find a mixed composition, containing a light mass fraction of ∼ 80%. Unless the extragalactic component becomes significant already below 10 17.5 eV, our measurements indicate an additional Galactic component dominating at this energy range.

Observations were made with the Low Frequency Array (LOFAR [START_REF] Van Haarlem | LOFAR: The LOw-Frequency ARray[END_REF] ), a radio telescope consisting of thousands of crossed dipoles, with built-in air shower detection capability [START_REF] Schellart | Detecting cosmic rays with the LOFAR radio telescope[END_REF] . LOFAR records the radio signals from air showers continuously while running astronomical observations simultaneously. It comprises a scintillator array (LORA), that triggers the readout of buffers, stor-ing the full waveforms received by all antennas.

We have selected air showers from the period June 2011 -January 2015 with radio pulses in at least 192 antennas. The total uptime was ∼150 days, limited by construction and commissioning of the telescope. Showers that occurred within an hour from lightning activity, or have a polarisation pattern that is indicative of influences from atmospheric electric fields are excluded from the sample [START_REF] Schellart | Probing Atmospheric Electric Fields in Thunderstorms through Radio Emission from Cosmic-Ray-Induced Air Showers[END_REF] .

Radio intensity patterns from air showers are asymmetric due to the interference between geomagnetic and charge excess radiation. They can be reproduced from first principles by summing the radio contributions of all electrons and positrons in the shower. We use the radio simulation code CoREAS [START_REF] Huege | Simulating radio emission from air showers with CoREAS[END_REF] , a plug-in of CORSIKA [START_REF] Heck | CORSIKA: a Monte Carlo code to simulate extensive air showers[END_REF] , which follows this approach.

It has been shown that X max can be accurately reconstructed from densely sampled radio measurements [START_REF] Buitink | Method for high precision reconstruction of air shower X max using two-dimensional radio intensity profiles[END_REF] . We use a hybrid approach, simultaneously fitting the radio and particle data. The radio component is very sensitive to X max , while the particle component is used for the energy measurement.

The fit contains four free parameters: the shower core position (x, y), and scaling factors for the particle density f p and the radio power f r . If f p deviates significantly from unity, the reconstructed energy does not match the simulation and a new set of simulations is produced. This procedure is repeated until the energies agree within uncertainties. The ratio between f r and f p should be the same for all showers and is used to derive the energy resolution of 32% (see Figure 1).

The radio intensity fits have reduced χ 2 -values ranging from 0.9 to 2.9. All features in the data are well reproduced by the simulation (see Extended Data Figs. 1-5), demonstrating that the radiation mechanism is now well-understood. The reduced χ 2 -values exceeding unity may indicate remaining uncertainties in the antenna response, atmospheric properties, or limitations of the simulation software.

Radio detection becomes more efficient for higher-altitude showers that have larger footprints. The particle trigger, however, becomes less efficient since the number of particles reaching the ground decreases. To avoid a bias, we require that all the simulations produced for a shower pass the trigger criteria. Above 10 17 eV this cut removes 4 showers from the sample. At lower energies, this number rapidly increases, and we exclude all showers below 10 17 eV from this analysis.

Furthermore, we evaluate the reconstructed core positions of all simulated showers. Showers with a mean reconstruction error above 5 m are rejected. This cut does not introduce a composition bias, because it is based on the sets of simulated showers, and not on the data. The final event sample contains 118 showers.

The uncertainty on X max is determined independently for all showers [START_REF] Buitink | Method for high precision reconstruction of air shower X max using two-dimensional radio intensity profiles[END_REF] , and has a mean value of 16 g/cm 2 (see Extended Data Figure 6). Figure 2 shows our measurements of the average X max , which are consistent with earlier experiments using different methods, within statistical un-certainties. The high resolution for X max per shower allows us to derive more information on the composition of cosmic rays, by studying the complete shape of the X max distribution. For each shower, we calculate:

a = X proton -X shower X proton -X iron (1)
where X shower is the reconstructed X max , and X proton and X iron are mean values predicted by the hadronic interaction code QGSJETII.04 [START_REF] Ostapchenko | QGSJET-II: results for extensive air showers[END_REF] .

The cumulative probability density function (CDF) for all showers is plotted in Fig. 3. First, we fit a 2-component model of proton and iron nuclei, with the mixing ratio as the only free parameter. To calculate the corresponding CDFs we use a parametrisation of the X max distribution fitted to simulations based on QGSJETII.04. The best fit is found for a proton fraction of 62%, but it describes the data poorly with a p-value of 1.1 × 10 -6 .

A better fit is achieved with a four-component model (p+He+N+Fe), yielding a p-value of 0.17. While the best fit is found for a Helium fraction of 80%, the fit quality deteriorates only slowly when replacing helium by protons. This is demonstrated in Figure 4 where the p-value is plotted for four-component fits where the fractions of helium and proton are fixed, and the ratio between N and Fe is left as the only free parameter. The total fraction of light elements (p+He) is in the range [0.38,0.98] at 99% confidence level, with a best fit value of 0.8. The heaviest composition that is allowed within systematic uncertainties still has a best fit p+He fraction of 0.6, and a 99% confidence level interval of [0.18, 0.82]. The online method section contains information about the systematic uncertainties and the statistical analysis.

The abundances of individual elements depend on the hadronic interaction model. The X max values predicted by EPOS-LHC [START_REF] Pierog | EPOS Model and Ultra High Energy Cosmic Rays[END_REF] are on average 15-20 g/cm 2 higher than QGSJETII.04 (see Fig. 2). This coincides with the separation between, for example, protons and deuterium or between helium and beryllium. We therefore prefer to present our result as a total fraction of light elements, instead of placing too much emphasis on individual elements.

Recent results for the Pierre Auger Observatory indicate that the cosmic ray composition at 10 18 eV, just below the ankle, can be fitted with a mixture of protons and either helium (QGSJET.II04)

or nitrogen (EPOS-LHC) [START_REF] Aab | Depth of maximum of air-shower profiles at the Pierre Auger Observatory. II. Composition implications[END_REF] . With decreasing energy, their proton fraction drops, while their helium (or nitrogen) fraction rises, down to the threshold energy of 7 • 10 17 eV. An extrapolation of this trend to our mean energy of 3 • 10 17 eV connects smoothly to our best fitting solution in which helium dominates.

KASCADE-Grande has reported an ankle-like feature at 10 17.1 eV, where the spectral index for light elements changes to γ = -2.79±0.08 4 . However, they find a light particle (p+He) fraction below 30% at 3 × 10 17 eV (based on their Figure 4), considerably lower than our value. In contrast to LOFAR and Auger, their composition measurements are based on the muon/electron ratio.

Auger has reported a muon excess compared to all commonly used hadronic interaction models [START_REF] Aab | Muons in air showers at the Pierre Auger Observatory: Mean number in highly inclined events[END_REF] .

Inaccurate predictions of muon production, or X max , can be the cause of the discrepancy in the fraction of light particles between LOFAR and KASCADE-Grande.

If the knee in the all-particle spectrum near 3 × 10 15 eV corresponds to the proton or helium cut-off of the main Galactic cosmic-ray population, the corresponding iron cut-off would lie at most at an energy 26 times larger. If this population still dominates at 10 17 eV, the mass composition should be dominated by heavy elements at that energy. Therefore, the large component of light elements observed with LOFAR must have another origin.

In principle, it is possible that we observe an extragalactic component. In that case the ankle in the cosmic-ray spectrum, slightly above 10 18 eV, does not indicate the transition from Galactic to extragalactic origin. Instead, it can be explained as the imprint of pair production on the cosmic microwave background on an extragalactic proton spectrum [START_REF] Aloisio | A dip in the UHECR spectrum and the transition from galactic to extragalactic cosmic rays[END_REF] . However, since this feature only appears for a proton-dominated flux it is in tension with our data that favours a mixture of light elements.

A second Galactic component, dominating around 10 17 eV, can be produced by a class of extremely energetic sources (Galactic exatrons), like the explosions of Wolf Rayet stars into their stellar winds [START_REF] Stanev | Cosmic rays. IV. The spectrum and chemical composition above 10 GeV[END_REF] , or past Galactic gamma-ray bursts [START_REF] Calvez | Role of Galactic sources and magnetic fields in forming the observed energy-dependent composition of ultrahigh-energy cosmic rays[END_REF] . Alternatively, the original Galactic population could be reaccelerated by the Galactic wind termination shock [START_REF] Jokipii | Ultra-high-energy cosmic rays in a galactic wind and its termination shock[END_REF] . Such scenarios predict mixtures of light elements, consistent with our results. The distribution of f r / f p is fitted with a Gaussian, yielding σ = 0.12 on a logarithmic scale, corresponding to an energy resolution of 32%. This value is actually the quadratic sum of the energy resolution of the radio and particle resolutions. In this analysis, there was no absolute calibration for the received radio power yet, so f r has an arbitrary scale. solving for the nitrogen/iron ratio. The p-value (see Fig. 3) is plotted as a function of the proton and helium fraction. The optimal fit is found at 0% proton and 79% helium (p=0.17 No. of simulations Extended Data Figure 7 | Energy reconstruction Distributions of the ratio between true and reconstructed energy for proton (blue, solid) and iron showers (red, dashed). The two types of showers have a systematic offset of the order of ∼ 1%.

Methods

Event selection Cosmic ray detection at LOFAR continuously runs in the background during astronomical observations. When 16 out of the 20 scintillator stations of the LORA particle array detect a signal, a trigger is issued and the ring buffers of all active antennas within a ∼ 1 km radius are stored for offline analysis [START_REF] Schellart | Detecting cosmic rays with the LOFAR radio telescope[END_REF] . Which antennas are active depends on the settings of the astronomical observation. For this analysis we have selected showers that were measured with at least 4 antenna stations (corresponding to at least 192 antennas) in the low band (30 -80 MHz after filtering)

The trigger and selection criteria introduce a composition bias. This bias is removed with a cut based on dedicated sets of simulations that are produced for each observed shower. These sets contain 50 proton and 25 iron showers, that span the whole range of possible shower depths. A

shower is only accepted when all simulations in its set pass the triggering and selection criteria.

This anti-bias cut removes many showers below 10 17 eV, but only 4 above that energy. In this analysis, we restrict ourselves to the higher-energy showers and impose a cut on energy E reco > 10 17 eV.

The energy cut itself is another potential source of compositional bias, as the reconstructed energy might be dependent on the depth of the shower. However, in our reconstruction approach this effect is very small because energy and X max are fitted simultaneously. Extended Data Fig-

ure 7 shows distributions of the ratio between true and reconstructed energy for proton and iron simulations. The systematic offset between the two particle types is of the order of ∼ 1%.

We used data from the Royal Netherlands Meteorological Institute to check for lightning storm conditions during our observations. When lightning strikes have been detected in the North of the Netherlands within an hour from a detection, the event is flagged and excluded from the analysis. The presence of electric fields in the clouds can severely alter the radio emission even in the absence of lightning discharges [START_REF] Buitink | Amplified radio emission from cosmic ray air showers in thunderstorms[END_REF] . The polarisation angle of the radio pulse is very sensitive to the nature of the emission mechanism [START_REF] Schellart | Probing Atmospheric Electric Fields in Thunderstorms through Radio Emission from Cosmic-Ray-Induced Air Showers[END_REF][START_REF] Schellart | Polarized radio emission from extensive air showers measured with LOFAR[END_REF] and is used as an additional veto against strong field conditions.

Finally, a quality cut is imposed on the sample in order to only include showers that have a core position and arrival direction that allows accurate reconstruction. We use the dedicated sets of simulations produced for each shower to derive uncertainties on core position, energy and X max .

These three values are highly correlated, so a single cut on the core uncertainty of σ core < 5 m is sufficient.

The quality cut is based on the dedicated sets of simulations. These sets are produced for a specific combination of core position and arrival direction. Therefore, the quality cut is effectively a cut on position and direction, and does not introduce a composition bias.

Furthermore, we stress that there is no cut on the quality of the reconstruction of the actual data. By applying the cuts described above we obtain a sample of 118 showers that are fitted to the simulation yielding reduced χ 2 -values in the range 0.9-2.9. Deviations from unity can be ascribed to uncertainties in antenna response, atmospheric properties like the index of refraction, or limitations of the simulation software.

Reconstruction The energy and X max of the shower are reconstructed with the technique described in Buitink et al. ( 2014) [START_REF] Buitink | Method for high precision reconstruction of air shower X max using two-dimensional radio intensity profiles[END_REF] .

Statistical uncertainty

The statistical uncertainty on the power measurements of individual antennas include three contributions. First, there is contribution from the background noise which is a combination of system noise and the Galactic background. Secondly, there is a contribution from uncertainties in the antenna response model. There can be differences between the responses of antennas, either because of antenna properties (e.g. cross-talk between nearby antennas), or because of signal properties (e.g. polarisation). Since these fluctuations are different for each shower core position and arrival direction, they are essentially random and included as a 10% statistical uncertainty on the power. A third contribution is due to the error introduced by interpolating the simulated pulse power. Strictly speaking this is not a measurement uncertainty, but it must be taken into account when fitting the data to simulation. The interpolation error is of the order of 2.5% of the maximum power [START_REF] Buitink | Method for high precision reconstruction of air shower X max using two-dimensional radio intensity profiles[END_REF] . The three contributions are added in quadrature and produce the one sigma error bars shown in Extended Data Figures 12345.

The statistical uncertainty on X max is given by the quadratic sum of the uncertainties due to reconstruction technique and the atmospheric correction. The former is found by applying our analysis to simulated events with added Gaussian noise, where the noise level is determined from the data.

In the CORSIKA simulations the standard US atmosphere model was used. The reconstructed shower depth is corrected for variations in the atmosphere using data from the Global Data Assimilation System (GDAS) of the NOAA National Climatic Data Center. We follow a procedure developed by the Pierre Auger collaboration [START_REF] Abreu | Description of atmospheric conditions at the Pierre Auger Observatory using the Global Data Assimilation System (GDAS)[END_REF] . This typically leads to adjustments of the order of 5-20 g/cm 2 . Remaining uncertainty after correction if of the order of 1 g/cm [START_REF] Murase | Testing the hadronuclear origin of PeV neutrinos observed with IceCube[END_REF] The index of refraction of air is a function of temperature, air pressure, and relative humidity.

Using local weather information the final data sample was split in two equal size groups corresponding to conditions with relatively high or low index of refraction. The mean reconstructed X max of these groups deviate from the mean by ±5 g/cm 2 , and we adopt this value as an additional statistical uncertainty. Because the refractivity used in simulation corresponds to dry air there is also an associated systematic error (see below).

The total statistical uncertainty on X max is found by adding above factors in quadrature. A distribution of the uncertainty for the showers in our final sample is shown in Extended Data Figure 6.

The energy resolution is 32% and is found by comparing energy scaling factors of the radio power and particle density fit (see Figure 1).

Systematic effects

The data has been subjected to several tests to find the systematic uncertainty on the reconstructed values for X max :

• Zenith angle dependence The final data sample is split into two groups of equal size by selecting showers with a zenith angle below or above 32 degrees. For both groups the mean reconstructed X max is calculated, yielding deviations from the mean value of ±8 g/cm 2 . This spread is larger than expected from random fluctuations alone and is included as a systematic uncertainty. The dependence on zenith angle may be related with atmospheric uncertainties (see below).

• Index of refraction of air As explained above, the index of refraction changes because of differences in atmospheric conditions. Fluctuations on X max due to changing humidity are of the order of 5 g/cm 2 with respect to the mean. However, the index of refraction that was used in the radio simulations corresponds to dry air, and is a lower bound to the actual value.

Therefore, the real value of X max can be higher than the reconstructed value but not lower, and we adopt an asymmetrical systematic uncertainty of +10 g/cm 2 .

• Hadronic interaction model Since the reconstruction technique is based on full Monte Carlo simulations, it is sensitive to the choice of hadronic interaction model that is used. It has been shown with a comparison between QGSJETII.04, SYBILL 2.1, and EPOS-LHC, that the uncertainty due to model dependence is ∼ 5 g/cm [START_REF] Murase | Testing the hadronuclear origin of PeV neutrinos observed with IceCube[END_REF] . Note that the uncertainty on the composition due to different models (in other words: on how to interpret the measured X max values) is of course larger.

• Radiation code For this analysis we have used the radiation code CoREAS in which the contributions of all individual charges to radiation field are added together. The advantage of this microscopic approach is that it is completely model-independent and based on first principles. ZHAireS 34 is another microscopic code and gives very similar results [START_REF] Huege | The convergence of EAS radio emission models and a detailed comparison of REAS3 and MGMR simulations[END_REF] . To calculate the emission CoREAS uses the end-point formalism [START_REF] James | General description of electromagnetic radiation processes based on instantaneous charge acceleration in "endpoints[END_REF] , while ZHAireS is based on the ZHS algorithm [START_REF] Zas | Electromagnetic pulses from high-energy showers: Implications for neutrino detection[END_REF] . Both formalisms are derived directly from Maxwell's equations and have been shown to be equivalent [START_REF] Belov | Accelerator measurements of magnetically-induced radio emission from particle cascades with applications to cosmic-ray air showers[END_REF] . The other difference between CoREAS and ZHAires is that they take the particle distribution from different air shower propagation codes (CORSIKA and AIRES respectively) that internally use different hadronic interaction models. Since the radiation formalisms themselves are equivalent, small differences between CoREAS and ZHAireS are most likely due to differences in the hadronic interaction models used to simulate the particle interactions. The choice of radiation code does therefore not introduce an additional systematic uncertainty on top of the uncertainty to hadronic interaction models that is already included. A comparison study with LOFAR data did also not show any evidence for a systematic offset between the codes and will be published in an upcoming paper.

The remaining small dependence of X max on zenith angle is possibly related to the index of refraction. Showers with different inclination angles have their shower maximum at different altitudes, and therefore different local air pressure and index of refraction. Therefore, increasing the index of refraction used in simulations will result in a zenith-dependent change in reconstructed X max .

This possibly removes the observed dependence of the composition on zenith angle. Correctly taking into account a complete atmospheric model for the profile of the refractivity of air is subject of further study. Here, we treat the effect conservatively by adding the first two contributions to the uncertainty linearly. The other two contribution are independent and are added in quadrature, yielding a total systematic uncertainty of +14/ -10 g/cm 2 .

The systematic uncertainty on the energy reconstruction with the LORA particle detector array is 27%, which includes effects due to detector calibration, hadronic interaction models, and the assumed slope of the primary cosmic-ray spectrum in the CORSIKA simulations [START_REF] Thoudam | LORA: A scintillator array for LOFAR to measure extensive air showers[END_REF][START_REF] Thoudam | Measurement of the cosmic-ray energy spectrum above 10 16 eV with the LOFAR Radboud Air Shower Array[END_REF] Statistical analysis For each observed shower, we calculate:

a = X proton -X shower X proton -X iron (2)
where X shower is the reconstructed X max , and X proton and X iron are mean values predicted by QGSJETII.04 [START_REF] Ostapchenko | QGSJET-II: results for extensive air showers[END_REF] . Thus a is an energy-independent parameter that is mass sensitive. A pure proton composition would give a wide distribution of a centered around zero, while a pure iron composition gives a narrower distribution around unity.

From the measurements we construct a cumulative distribution function (CDF) in the following Monte Carlo approach. A realisation of the data is made by taking the measured values for the energy and X max , adding random fluctuations based on the statistical uncertainty of these parameters, and calculating the a parameters and the corresponding CDF. By constructing a large number of realisations with different random fluctuation, we can calculate the mean CDF and the region that contains 99% of all realisations. These are indicated in Figure 3 as the solid blue line and the shaded region respectively.

We fit theoretical CDFs based on composition with two or four mass components to the data.

The test statistic in the fit is the maximum deviation between the data and the model CDFs. The p-value is given by the probability of observing this deviation, or a larger one, assuming the fitted

Figure 1 |

 1 Figure 1 | Energy resolution. The distribution of f r / f p is fitted with a Gaussian, yielding σ = 0.12

Figure 2 |Figure 3 |Figure 4

 234 Figure 2 | Measurements of X max . The mean depth of shower maximum as a function of en-

Extended Data Figure 1 |Extended Data Figure 2 |Extended Data Figure 3 |Extended Data Figure 5 |Extended Data Figure 6 |

 12356 ), but the deviation only deteriorates slowly when replacing helium with proton. The thick black contour line contains all combinations for which p > 0.01. At this significance level the total fraction of light elements (p+He) lies between 0.38 and 0.98. Fitted lateral distributions. Lateral distribution of radio pulse power for all 118 measured showers (red circles) and the corresponding best-fitting CoREAS simulation (blue squares). Fitted lateral distributions. Continuation of Extended Data Figure1. Fitted lateral distributions. Continuation of Extended Data Figure2. Fitted lateral distributions. Continuation of Extended Data Figure4. Distribution of uncertainty on X max The distribution of the uncertainty on X max for all showers used in this analysis. The mean value is 16 g/cm 2 .

Acknowledgements We acknowledge financial support from the Netherlands Organization for Scientific Research (NWO), VENI grant 639-041-130, the Netherlands Research School for Astronomy (NOVA), the Samenwerkingsverband Noord-Nederland (SNN) and the Foundation for Fundamental Research on Matter (FOM). We acknowledge funding from the European Research Council under the European Union's Seventh Framework Program (FP/2007-2013) / ERC (grant agreement n. 227610) and under the European Union's Horizon 2020 research and innovation programme (grant agreement n. 640130). LOFAR, the Low Frequency Array designed and constructed by ASTRON, has facilities in several countries, that are owned by various parties (each with their own funding sources), and that are collectively operated by the International LOFAR Telescope (ILT) foundation under a joint scientific policy.

Author Contributions All authors are part of the LOFAR collaboration and have contributed to the design, construction, calibration, and maintenance of LOFAR and/or LORA. The first thirteen authors constitute the Cosmic Ray Key Science Project and have contributed to the acquisition, calibration, and analysis of cosmic ray radio data and LORA data. The manuscript was written by S.B. and subjected to an internal collaboration-wide review process. All authors approved the final version of the manuscript.

Author Information Reprints and permissions information is available at www.nature.com/reprints. The authors declare no competing financial interests. Correspondence and requests for materials should be addressed to Stijn Buitink (Stijn.Buitink@vub.ac.be). composition model.

We first use a two-component model of proton and iron nuclei, where the mixing ratio is the only free parameter. The best fit is found for a proton fraction of 62%, but it describes the data poorly with a p-value of 1.1 × 10 -6 .

A better fit is achieved with a four-component model (p+He+N+Fe), yielding a p-value of 0.17. While the best fit is found for a Helium fraction of 80%, the fit quality deteriorates only slowly when replacing helium by protons. This is demonstrated in Figure 4 where the p-value is plotted for four-component fits where the fractions of helium and proton are fixed, and the ratio between N and Fe is left as the only free parameter. The solid line in this Figure contains the parameter space where p > 0.01. We construct a 99% confidence level interval on the total fraction of light elements (p+He) by finding the two extreme values of this fraction that still lie within the p > 0.01 region.

The total fraction of light elements (p+He) is in the range [0.38,0.98] at 99% confidence level, with a best fit value of 0.8. The heaviest composition that is allowed within systematic uncertainties (see above) still has a best fit p+He fraction of 0.6, and a 99% confidence level interval of [0.18, 0.82].