Dr Alexis A Boileau

Aïmane Cheikh

Arnaud Fouchet

Adrian David

Dr Christophe Labbé

Dr Philippe Marie

Dr Fabrice Gourbilleau

Ulrike Luders
email: ulrike.luders@ensicaen.fr

Dr A Cheikh

Dr A Fouchet

Dr U David

Crismat Lüders

Tuning of the Optical Properties of the Transparent Conducting Oxide SrVO3 by Electronic Correlations

Keywords: correlated materials, epitaxy, thin films, transparent conducting oxides, vanadates

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

Tuning of the Optical Properties of the Transparent Conducting Oxide SrVO 3 by Electronic Correlations

Alexis Boileau, Aïmane Cheikh, Arnaud Fouchet, Adrian David, Christophe Labbé, Philippe Marie, Fabrice Gourbilleau, and Ulrike Lüders* of SVO as a TCO is the fact that although the charge density is that of typical metals (in the range of 10 22 cm -3 ), the screened plasma frequency ω p is pushed out of the visible range due to the strong electronic correlations enhancing the effective mass m*. Consequently, the electronic correlations do not only govern the electronic transport properties, but also the optical ones in the visible range.

Such interplay between the electronic and optical properties via the electronic correlations may also open a new way for the optimization of the functional properties of SVO. In strongly correlated metals, it is known that strain may influence the correlation characteristics via the related crystal distortions, [3][4][5] which were also predicted theoretically for SVO. [6,7] In a strongly correlated TCO as SVO, it may therefore be possible to not only influence the electrical transport by strain, but also the optical properties. Although the structural and electronic properties of SVO were studied intensely in the thin film [8][9][10][11][12] or in superlattices form, [10,13] no comparative study of SVO films on different substrates inducing variation of strain was published, neither its influence on the optical properties. The aim of this study is therefore to use three substrates, on which epitaxial growth of SVO is possible, in order to elucidate the influence of the strain on both the electronic and the optical properties.

However, in the discussion of the relationship between structure and properties in transition metal oxide thin films, the stoichiometry of the films has also to be taken into account. Especially in SVO, it has been shown that both a cationic [9] as also an anionic nonstoichiometry [11,12,14,15] play an important role in the properties. The management of the oxygen content is a special issue in SVO due to the existence of some oxides with higher oxidation state of vanadium such as Sr 3 V 2 O 8 or Sr 2 V 2 O 7 , being thermodynamically more stable than the perovskite structure. [16] This implies that widely used mechanisms to influence on the oxygen stoichiometry of thin films, as the growth in different oxygen pressures or postannealing treatments, are not possible or have to be controlled closely in order to avoid the formation of the orthovanadate, which is transparent but insulating. Even more, recently it has been shown in a closely related system, SrCoO 3 , that strain and oxygen vacancies show a complex interplay: [17] in films undergoing a tensile strain, the density of oxygen vacancies is enhanced,

The vanadate SrVO 3 is a transparent conductor perovskite with optical and electrical properties competing with those of the most-used indium tin oxide material. Although its charge density is comparable to that of metals, SrVO 3 shows a plasma frequency below the visible range due to strong electronic correlations characterizing the electronic transport in this material and enhancing the effective mass. Therefore, the well-known interplay between the structure and the electronic properties of strongly correlated systems can be used in such transparent conductor to tune the optical properties, as the plasma frequency also depends on the effective mass. In this study, SrVO 3 films are grown by pulsed laser deposition onto different lattice mismatched perovskite substrates such as SrTiO 3 , LaAlO 3 , and (LaAlO 3 ) 0.3 (Sr 2 TaAlO 6 ) 0.7 at different growth temperatures. The structural, electronic, and optical properties are analyzed, illustrating the influence of the strain on the structure of the films and on a shift of the plasma frequency. The electronic correlations in this new group of transparent conducting oxides can be therefore used as a supplementary lever for the tuning of the functional properties.

Introduction

SrVO 3 (SVO) is one of the strongly correlated metals that were recently identified to be transparent conducting oxides (TCO). [1] The demonstrated properties of this system stand up to those of the standard indium tin oxide material, providing therefore a valuable alternative for a number of applications. Although it was shown that the amorphous phase of SVO is not conducting, [2] the potential of application stays important in adapted systems due to the absence of indium and its price fluctuations related to the limited availability. An intriguing aspect probably due to the weakening of the CoO bonds. This observation seems to be part of more general tendency in perovskite oxides, where the introduction of oxygen vacancies is favored by tensile strain allowing for an accommodation of the oxygen vacancy volume. [18] In order to estimate the effects of oxygen vacancies, which cannot be measured with easy means in thin films, we have included also a study on the growth temperature T G . The films, being deposited in vacuum to avoid in situ overoxidation, undergo a strongly reducing atmosphere at least at high T G , while at low T G the oxygen transport in the film is less activated. By comparing the films at different T G , it will be possible to elucidate the role of oxygen vacancies on the functional properties of the film, although it will not be possible to gain absolute knowledge.

In this work, we report on the experimental investigation of the growth of SrVO 3 films on different perovskite substrates and their influence on the structural, electrical, and optical properties. The selected substrates are (001) oriented SrTiO 3 (STO), (LaAlO 3 ) 0.3 (Sr 2 TaAlO 6 ) 0.7 (LSAT), and LaAlO 3 (LAO) monocrystals, adapted to the growth of SVO. The lattice parameter of STO, LSAT, and LAO substrates is 3.905, 3.868, and 3.790 Å, which lead to a lattice mismatch with SVO bulk (3.840 Å) of 1.69%, 0.73%, and -1.30%, respectively. Therefore, SVO films grown on STO and LSAT are tensile-strained while the films grown on LAO are compressive-strained. The aim of this study is to determine the functional properties for films undergoing the different strain states and understand the role of these strains.

Results and Discussion

Structural Characterizations

First, the picture of deposited films on STO, LAO, and LSAT displayed in Figure 1a gives an overview of the good transparency of SVO films for a thickness of 30 nm. Concerning X-ray diffraction (XRD) measurements, Figure 1b shows specular θ-2θ diffractograms around the (002) Bragg peak of SVO films grown on STO and LAO between 300 and 700 °C (the data concerning the films on LSAT were published elsewhere [2] ). At 300 °C, only the peak attributed to the (002) reflection arising from the substrate is visible at 46.47° and 47.97° for STO and LAO, respectively. The peak due to the (002) reflection of the SVO cubic structure appears around 47.25°, and 46.50° on STO and LAO, respectively, between 300 and 400 °C. Similar results were found for the films on LSAT. This crystallization temperature of SVO is remarkably low and does not seem to depend on the substrate, at least in the defined temperature interval. From the observed Laue fringes, the layer thickness was determined to be in the range of 37-40 nm by simulations (not shown). As the presence of Laue fringes can be considered as a marker of the quality of the crystalline local order, highly crystallized films are grown at 400 °C on STO and LSAT with a flat film/substrate interface. By contrast, the film deposited on LAO at the same temperature does not exhibit Laue fringes. In the case of LSAT, the lattice mismatch is limited to 0.72% but increases positively up to 1.66% for STO and negatively down to -1.32% for LAO. Given that the epitaxial tensile strain induced by the STO substrate is more important than the compressive one induced by LAO, the loss of oscillations may indicate the negative impact of compressive strains on the crystallization of SVO. The energy required for the epitaxial stabilization of SVO in compressive stress is therefore partially reinvested in the elastic potential energy required to induce the deformation of the perovskite structure by the apical elongation of the octahedra of the corner-sharing perovskite structure.

As T G is increased, the amplitude of the Laue fringes around the (002) SVO peak deposited on LAO significantly increases as the provided thermal energy becomes high enough to ensure both the deformation and the crystallization of the structure. Consequently, the quality of the crystallization of SVO at low temperature can be strongly modulated by the deformation of the strained epitaxial structure imposed by

(b)

Figure 1. a) Image of SVO films grown onto STO, LAO, and LSAT between 300 and 700 °C. b) θ-2θ X-ray diffractograms in the vicinity of the (002) reflection of SVO thin films on STO and LAO grown between 300 and 700 °C. For SVO/LSAT samples refer to ref. [2].

the substrate and widely enhanced by implication of a biaxial tensile strain.

To get better insight in the in-plane structural relation between SVO and the substrates, asymmetric film reflections were also measured in skew symmetry. In Figure 2a, reciprocal space maps (RSMs) around the (103) Bragg peaks for the film on STO are presented for T G at 400 and 700 °C, representative also for the other used substrates. From the vertical alignment of the film and substrate peaks at the same Q x value, we can note that SVO films are fully strained to the substrate regardless of growth temperature for STO and slightly relaxed for LAO.

From the (002) and (103) peaks, we have determined the evolution of in-plane and out-of-plane cubic cell parameters a f and c f of SVO films depending on the growth temperature. The corresponding lattice parameters are displayed in Figure 2b. The in-plane lattice parameter a f (open symbols) of SVO is found to be 3.90(9)-3.91(5), 3.86(2)-3.86(5), and 3.80(1)-3.79(5) Å at 400 and 700 °C on STO, LSAT, and LAO, respectively, whereas one of the substrates a s is equal to 3.90(5), 3.86 (8), and 3.79(0) Å for bulk substrates. Therefore, as observed earlier in the RSMs, a f matches closely the substrate lattice parameter for all films, independently of the induced strain. However, the out-of-plane lattice parameter c f (solid symbols) exhibits a different evolution with T G and depends strongly of the nature of the substrate. For SVO//STO, c f increases rapidly from 3.84(5) to 3.87(5) Å as the T G increases up to 600 °C. It can be noticed that for the film deposited at 400 °C, the c-axis is close to the value reported for bulk SVO, found to be between 3.82 and 3.84 Å. [19][20][21][22] Around 600 °C, c f goes through a maximum and decreases slightly down to 3.86 Å at 700 °C. The expansion of the c-axis observed between 400 and 600 °C may be related to the creation of oxygen vacancies in SVO due to the growth under vacuum. The incorporation of oxygen vacancies induces a swelling effect of the cell driven by the reduction of V 4+ to the larger cation V 3+ : the ionic radius increases slightly from 0.58 to 0.64 Å. [START_REF] Lide | CRC Handbook of Chemistry and Physics: A Ready-Reference Book of Chemical and Physical Data: 2006-2007[END_REF] Actually, the near-stagnation of the c-axis of SVO//STO at high T G (whereas it increases continuously for LAO and LSAT) can be understood by the vicinity of the STO substrate. Heating of the STO under vacuum is known to generate oxygen vacancies in the structure and to create colorcenters. [START_REF] Frye | [END_REF][25][26] Conversely, SVO is highly sensitive to oxygen capture even for reduced oxygen partial pressures of 10 -4 Pa. [12] At T G above 500-600 °C, the thermal energy is high enough to activate the oxygen diffusion at the interface from STO toward SVO and compensate the oxygen depletion in the film created by the reducing growth conditions. Consequently, the substrate can be viewed as an oxygen reservoir for the film during the growth. Thus, a change of T G would lead to a modification of the vacancy doping in STO and SVO, the film being a source of oxygen vacancies, as proposed also by Herranz et al. [27] For SVO//LAO films, c f is remarkably shifted toward high values to accommodate compressive strains induced by the LAO misfit. However, LAO is stable under vacuum conditions and the reduction of the substrate does not occur as much as in STO one. In these conditions, we observe that c f rises continuously from 3.904 to 3.925 Å when the deposition temperature is increased, indicating the creation of oxygen vacancies in SVO for higher T G . At the same time, the crystallinity of the films is improved due to the increased thermal energy.

For SVO//LSAT, we observe a similar evolution of c f with T G . The c-axis lattice parameter increases from 3.905 to 3.915 Å as observed on the other substrates, but the high observed values, similar to SVO deposited on LAO, is somewhat unexpected due to the change of the misfit from compressive on LAO to slightly tensile on LSAT. If the reaction of the film on the strain would be purely elastic, values smaller than the bulk would be expected on LSAT. It is known that a deviation in the cation stoichiometry of the perovskite could influence the lattice parameters of films. Studies of SVO thin film on LSAT substrate have also variations between 3.82 and 3.84 Å depending on the Sr and V content. [9] The reported values are much smaller than those observed in this study, and Rutherford backscattering spectrometry (RBS) showed a Sr/V ratio of about 1, so that this effect can be ruled out here. However, it was shown that the presence of biaxial tensile stain may induce an oxygen depletion to satisfy the deformation of the octahedral network. Recently, Petrie et al. [17] have shown a similar behavior in SrCoO 3 films, where the creation of oxygen vacancies was shown to be facilitated by the introduction of tensile strain, probably due to weakened MO bonding. Therefore, coherently strained films were observed on all substrates, permitting to study the influence of the strain on the physical properties. However, care has to be taken with the films on STO, where an important influence from possible oxygen vacancies cannot be excluded.

Transport Properties

Transport Properties of As-Grown Films

The temperature dependence of the electrical resistivity of the films was measured between 5 and 300 K. Results are reported in Figure 3 for STO, LAO and LSAT at different T G . For all crystallized films, we observe a metallic behavior in the whole temperature range in good agreement with earlier works. [1,15] Globally, the resistivity measured at room temperature is between 1.10 and 2.37 × 10 -4 Ω cm and corresponds to the typical order of magnitude observed for SVO films deposited by pulsed laser deposition (PLD). [START_REF] Ritums | AIP Conf. Proc[END_REF][START_REF] Yip | [END_REF] Despite large differences observed in the structural quality of films depending on the growth temperature and the nature of the substrate, the resistivity remains in the same order of magnitude and does not vary much at room temperature. At most, we observe a varia-Figure 3. dependence of the electrical resistivity of SVO films grown on STO, LAO, and LSAT substrates between 400 and 700 °C. Data in the bottom panel are reproduced from Boileau et al. [2] tion in the resistivity of a factor of two. Hence, the influences of the substrate and T G are relatively small compared to the need of well-defined crystalline perovskite structure for good transport properties, as was also pointed out for thin films on LSAT. [2] However, regarding the details of the transport behavior, differences due to the substrates can be observed. In Figure 4, different markers for all substrates and the different growth temperatures are combined. In Figure 4a, the electrical resistivity measured 5 K (ρ 5K ) as well as the residual resistivity ratio (RRR) defined as the ratio between ρ 300K and ρ 5K are displayed. The temperature dependence of the resistivity was fitted on the basis of the typical model for metallic conductors ρ(T) = ρ 0 + A × T α where ρ 0 is the residual resistivity taken at 0 K and A is related to the electronic interactions for a strongly correlated metal. The exponent α gives an indication on the Fermi liquid character of the system. [10,15,30] The fit results are shown in Figure 4b. To complete electrical transport data, the density and the mobility of carriers measured by Hall effect at 300 and 5 K are reported in Figure 4c,d,respectively.

For the films on LAO and LSAT, the evolution of the different markers with T G is remarkably similar. The resistivity decreases in both cases with the increase of T G , leading to an increase of the RRR. For high T G , the fit parameters of the temperature dependence of the resistivity are near to what was observed in the bulk material, [31] α = 2 and A = 5.9 × 10 -10 Ω cm K -2 . Regarding the characteristics of the charge carriers, the charge density N decreases with T G , reaching a value about between 1.5 (5 and 1.7 × 10 22 cm -3 (300 K) for the samples grown at 700 °C. This value is consistent with the expected value based on the hypothesis that each V 4+ ion (3d 1 ) contributes one electron, which would add up to N = 1.77 × 10 22 cm -3 with the bulk density. The mobility µ increases with T G , reaching values of 6-7 cm 2 V -1 s -1 at 5 K for T G = 700 °C. These values are in the same range as those observed in other studies of thin films on LSAT substrates, outreaching those of other films grown by PLD, [10] but staying about a factor of two (at room temperature) lower than those of molecular beam epitaxy (MBE) grown samples. [1,9] The films on LAO and LSAT show therefore charge transport characteristics, which are comparable to bulk properties at high T G , while for low T G , the markers related to disorder (ρ 5K , RRR, α, and µ) show stronger deviations. Thus, the observed evolution of the resistivity with T G can be explained by the decrease of disorder and therefore to the increase in crystalline quality, as was observed in the XRD study on LAO. However, for the films grown on LSAT, an exceptionally high crystalline quality is already observed at 400 °C, [2] comparable to the one of the film grown at 700 °C, which does not seem consistent with structural disorder being at the origin of the degradation of the electronic transport. Moyer et al. [9] have shown that an imbalance of the Sr-V stoichiometry may also be at the origin of high residual resistivities and low mobilities, but the RBS analysis of the films of this study did not show an evolution for different T G or different substrates. The origin of the observed scattering may therefore also be related to non-crystalline disorder, i.e., purely electronic disorder, which will be addressed in future studies of the system.

The films grown on STO show a different behavior compared to the films grown on LAO and LSAT. While for the sample grown at 400 °C, the resistivity is still comparable to the films grown in LAO and LSAT, ρ increases for higher T G and the RRR decreases. This evolution is also observed for A and α, where the characteristics are comparable to the films on LAO and LSAT at 400 °C, but while these latter approach the bulk properties at higher T G , the films on STO deviate more strongly from the bulk values, with α decreasing down to 1.77 at T G = 600 °C. A increases in the same temperature range, resulting in a value about four times larger than the bulk value. Interestingly, the film grown at 700 °C reverses the tendencies for the low T G samples, probably due to the starting exchange of the oxygen vacancies between the substrate and the film, as was observed in the structural properties. N stays roughly constant in the observed T G range, and µ decreases, indicating that it is not a lack of charge carriers, but rather a lack of mobility leading to the high resistivity at high T G .

It seems therefore that the strain plays only an indirect role for the electronic properties, as the films on LAO (large compressive strain) and LSAT (small tensile strain) show similar behaviors, while the properties of the films on STO (large tensile strain) are governed rather by the oxygen stoichiometry in the film. This interpretation is underlined by the fact that the samples grown at 400 °C show similar properties on all substrates, at a growth temperature where the thermodynamic reduction of the film during the growth should be rather small. If the strain would have a direct influence, these films, fully strained, should show pronounced differences. On the other hand, the strain seems to have an indirect influence via the oxygen vacancies: for increasing T G , the probability of the introduction of oxygen vacancies in the films becomes more important for all films, but the ones grown on LAO and LSAT show actually a tendency to more bulk-like properties, excluding therefore a high oxygen vacancy density. For STO, the introduction of oxygen vacancies is favored both by the high tensile strain and by T G up to 600 °C, while at higher temperature the oxygen exchange between the substrate and the film sets in. The transport characteristics reflect this evolution, illustrating the importance of the oxygen vacancies and therefore indirectly of the tensile strain for the films on STO.

Turning back to the initial motivation of the study, the influence of the substrate on the electronic correlations can be estimated using the fit parameter A and the charge density N. In a simplified model, [31] A 1/2 N -1/3 is proportional to the effective electronic mass m*, one of the governing parameters also for the optical properties. In Figure 5, the values of the coefficient A 1/2 N -1/3 for the different substrates and T G are shown, confirming the observed differences in the transport properties of the films grown on STO and those grown in LAO and LSAT. While at T G = 400 °C, the values are similar for all substrates, A 1/2 N -1/3 enhances for films grown on STO for higher T G , while they decrease for films grown on LAO and LSAT. The strongest difference is observed for the samples grown at T G = 500 °C, with a ratio in the order of 2, indicating a value of m* twice as high on STO compared to the other substrates. The observed structural changes have therefore an influence on the electronic correlations of the SVO films.

Optical Properties of As-Grown Films

Optical transmission spectra were acquired in the visible range. The transmittance taken in the middle of the visible window at 500 nm is shown in Figure 6a for the films grown on STO and LSAT. Unfortunately, the acquisition of transmission spectra of SVO//LAO was not possible due to the presence of twined domains in the LAO substrates and the difficulty to use a reliable reference between the pristine substrate and the deposited one to realize the background in the dual beam acquisition mode.

The analysis of the transmittance shows a very high optical transparency for SVO films deposited on both STO and LSAT, even outreaching those of earlier reported data on LSAT substrates by Zhang et al. [1] It is interesting to note here, that the films grown on STO show an important transparency at low T G of nearly 93% going down to 80% in films grown at 700 °C. This evolution may be related to the creation of active colorcenters in the STO substrate with the presence oxygen vacancies. [25,32] In fact, it is not possible to distinguish between a loss of transparency between the STO substrate and the SVO films deposited at high temperature, so that the here reported values can only be considered as minimum values for the SVO films.

Regarding the reflectivity spectra of the SVO films shown in Figure 6b, well-known structures related to the free charge carriers (at high wavelength) and to a SVO intraband transition within the V 3d conduction band located at around 3.5 eV [21] (at low wavelength) can be identified. To get more into detail, these two signatures allow to determine the screened plasma frequencies ω p of the films by the observed minimum. [33] For some of the samples, ω p was also determined by the fitting of ellipsometry data, showing the same values and confirming therefore that the minimum of the reflectivity spectra is a good measure of ω p . The obtained values of the films on the different substrates are shown in Figure 7. It can be clearly observed that the films on STO show higher ω p than the films on LSAT and LAO. There is also a certain offset between these two latter substrates, although much smaller compared to STO. Taking into account the estimation of the effective mass m* from the electrical measurements (Figure 5), this dependence of ω p on the substrate type can be therefore understood by the changes of the electrical properties of the SVO films. In principle, it would be also possible to determine m* from ω p in TCOs, [34,35] but the relative permittivity has to be taken into account, which may depend also on the details of the growth process. Therefore, we refrain here from this calculation. However, the variation reaching around 10% of ω p on different substrates is reasonable regarding the variation of around a factor of two for the estimated variation of m* based on the electronic properties: A is proportional to m* 2 , while ω p is proportional to m* -1/2 .

Soft-Annealing Experiments on STO

In the goal to highlight the effect of the oxygen vacancies on the properties, additional SVO films grown on STO under the same deposition conditions have been annealed under air during 20 h at 100 °C. The annealing temperature was chosen carefully to avoid the massive oxidation of SVO to Sr 3 V 2 O 8 occurring above 300 °C. [16] Figure 8a shows the evolution of the out-of-plane lattice parameter determined by high-resolution XRD measurements of SVO before and after annealing.

After soft annealing under ambient atmosphere, the out-ofplane lattice parameter of SVO//STO decreases to approach the theoretical value of bulk SVO expected at 3.84 Å (represented by dashed line in Figure 8a). The swelling effect driven by the vanadium reduction in the presence of oxygen vacancies is compensated by the partial reoxidation of depleted sites under quite low temperature. The absolute variation of c f reported in Figure 8b shows that the out-of-plane deformation of the lattice is more important for films deposited at high T G . Hence, the depletion in oxygen in SVO films grown on STO under vacuum is increased by the deposition temperature.

Electrical characterizations of films summarized in Table 1 for the two uttermost T G show a significant improvement of the transport parameters after thermal treatment. Indeed, the resistivity decreases with a variation of 0.09 and 0.55 × 10 -4 Ω cm at 400 and 700 °C while the residual resistivity is diminished of 0.05 and 0.55 × 10 -4 Ω cm. At the same time, the residual resistivity ratio is nearly stable for 400 °C sample and increases notably from 1.22 up to 1.30 at 700 °C. In other words, the SVO sample grown at 400 °C is less sensitive to the thermal treatment than the one deposited at 700 °C.

This discrepancy can be related to an easier filling of oxygen vacancies during the treatment and consequently, a higher level in the oxygen depletion in SVO deposited at high temperature. This result is coherent with the evolution observed previously in the structural and electrical properties of SVO//STO and more precisely, the diminution of RRR indicating the oxygen loss with the growth temperature. We can notice that the carrier density and the mobility increase slightly after annealing. This variation is also more important at 700 °C compared to 400 °C. Although the mobility remains less than values expected around 5-10 cm 2 V -1 s -1 for bulk SVO or films deposited by MBE, [1] its increase reveals the reinforcement of the conduction pathway required for the electron delocalization in the structure via the O2p-V3d exchange.

The variation of the A coefficient before and after annealing is reported in Figure 9a for these samples. The fit parameter A increases with the growth temperature, in good agreement as was previously observed for the first series of SVO//STO samples (Figure 4b). A is dependent on two contributions: the average scattering probability and the electronic effective mass. [36] As a consequence, the increase of A with T G confirms the higher structural disorder in SVO//STO films deposited at high temperature. At the same time, the effective mass increases with T G and can be related to reinforced electronic correlations at high temperature. After soft thermal treatment, we observe a spectacular depletion of A with T G whereas it remains quasi-constant with no significant variations at 400 °C. The effect of the thermal treatment on m* is investigated by the plot of A 1/2 N -1/3 reported in Figure 9b. At 400 °C, the estimation of m* is stable but is divided by a factor of four at 700 °C being comparable to those on LAO and LSAT estimated around 1.1 × 10 -15 SI (see Figure 5). Figure 10 reports the reflectivity spectra of SVO films grown on STO at 400 and 600 °C before and after soft annealing at 100 °C during 20 h. First, a depletion is observed in the reflectivity level of about 5% in the entire range of the probed wavelengths, i.e., including both the interband transitions and the Drude contribution. This effect, constant through the entire visible spectrum, can be related to a modification of the extreme surface of films with the possible creation of an overoxidized layer or nanostructures such as Sr 2 V 2 O 7 evidenced elsewhere by Berini et al. when SrVO 3 is annealed in oxidizing conditions. [12] Second, the annealing leads to a displacement of λ p (ω p *) from 522 to 527 nm (1.19 to 1.18 eV) for sample grown at 400 °C and from 525 to 541 nm (1.18 to 1.15 eV) for the 700 °C one. This tendency is not consistent with the tendency of A 1/2 N -1/3 reported in Figure 9b, showing a decrease of m* after annealing. However, it can be observed that the position of the interband transitions situated between 350 and 400 nm moves to high values, having a direct influence on ω p *. The increase of ω p * after annealing is therefore rather related to the effects on the lattice probed by the interband transitions, rather than the properties of the charge carriers above. Thus, conclusions about m* on the optical measurements are not evident in this case.

Conclusion

In conclusion, SrVO 3 thin films have been grown by PLD on STO, LSAT, and LAO perovskite substrates and their structural, electrical and optical properties have been examined. First, SrVO 3 crystallizes surprisingly from 400 °C on whatever the used substrate with however a remarkably high structural quality in the case of in-plane tensile strain while a compressive strain limits the local order in the crystalline structure.

Second, the resistivity of films remains around (≈1-2) × 10 -4 Ω cm, indicating that transport properties are not too dependent on the substrate and the growth temperature, even if slight variations are observed. In particular, we have shown the depletion of oxygen vacancies in SVO deposited onto STO substrate at high temperature by postannealing experiments. The sub-stoichiometry in oxygen can be understood by an interdiffusion effect at the SVO/STO interface during the deposition under reducing conditions where STO is unstable.

Globally, the properties of SVO are not affected strongly by strains induced by the substrate, oxygen depletion in films seems to play a more important role on the conduction. However, the changes in the electronic properties indicate a higher effective electronic mass m* for the films grown in STO, resulting also in a higher screened plasma frequency ω*, compared to the other substrates.

In the objective to design and control tailored optical properties ensuring a good transparency through the visible range, an adapted system based on SVO film has to show a low plasma frequency, out of the visible range. As the plasma frequency is proportional to (N/m*) 1/2 , this condition can be accomplished by increasing m* by introducing, for example, oxygen vacancies in the structure via a chemical exchange with the substrate or by the tensile strains. However, important oxygen depletion affects the electronic transport of films and a compromise between the conduction and the transparency should be found.

Experimental Section

SVO thin films have been grown on double side polished (001) oriented STO, LSAT, and LAO substrates by PLD. Prior to deposition, substrates (5 × 5 × 0.5 mm) were successively cleaned in acetone and ethanol ultrasonic baths. A KrF excimer laser (λ laser = 248 nm) was used with a repetition rate of 3 Hz and a laser fluence of ≈1.6 J cm -2 focusing on a Sr 2 V 2 O 7 polycrystalline target prepared by standard solid-state reaction. The deposition was carried out under vacuum to induce the reduction of the target material and to avoid overoxidation of the SVO films, known to lead to the presence of higher oxides as Sr 3 V 2 O 8 . [12] The residual deposition pressure at the growth temperature was 1 × 10 -6 mbar, while the base pressure is 1 × 10 -7 mbar. The target-substrate distance was set at 50 mm in order to achieve a deposition rate of 0.09 Å per pulse and the SrVO 3 thickness was fixed to be around 40 nm for all samples by adjusting the deposition time. In order to investigate the influence of the temperature on the crystallization of SrVO 3 films, the substrate temperature (T G ) has been varied from 300 to 700 °C.

The structural properties and the thickness of the samples were characterized by XRD with a Bruker D8 Discover diffractometer operating with monochromatic Cu K α1 radiation (λ = 1.5406 Å). The temperature dependence of transport properties was determined by the four probe method from 5 to 300 K using the Van der Pauw configuration in a physical properties measurement system by Quantum Design. To determine the carrier concentration and the mobility of carriers, Hall measurements were performed in the same temperature range by applying perpendicularly to the film plane a magnetic field within a range from -9 to +9 T. The optical properties were measured in the UV-vis-near infrared range using a Lambda 1050 Perkin-Elmer spectrophotometer operating in reflection mode by means of the universal reflectance accessory and transmission mode. For all substrates, a pristine substrate was taken as reference. In order to ensure the absence of a possible deviation in the cationic contents of SVO films, chemical analyses have been performed by RBS with 1.7 MeV He + ions (not shown here). The quantification of the Sr/V ratio indicates a maximum variation from the stoichiometric value Sr/V = 1% of 6% between 400 and 700 °C without a clear dependence on the growth temperature or the used substrate. It thus confirms an excellent transfer of the cation stoichiometry of the target under the deposition conditions used. This discards a possible effect of the cationic composition on the structural or functional properties of the films.

Figure 2 .

 2 Figure 2. a) Reciprocal space maps around the (103) cubic reflection of SVO thin films grown at 400 and 700 °C on STO and LAO substrates. The reflections marked by a star (400 °C) are artefacts due to the diffractometer setup. b) Out-of-plane (solid symbols) and in-plane (open symbols) lattice parameters evolution as a function of the growth temperature of SVO films on STO (red squares), LSAT (blue triangles), and LAO (red circles).

Figure 4 .

 4 Figure 4. Temperature dependence of transport properties of SVO film grown on STO (green squares), LAO (red circles), and LSAT (blue triangles). a) Electrical resistivity at 5 K and the RRR ratio, b) A and α coefficients, c) the density of free carriers at 300 and 5 K, and d) the mobility of free carriers at 300 and 5 K.

Figure A 1 / 2 N

 12 Figure A 1/2 N -1/3 proportional to m* as a function of growth temperatures T G for different substrates.

Figure 6 .

 6 Figure 6. a) Transmittance at 500 nm as a function of the substrate temperature and b) reflection spectra of SVO films in the vis-NIR range for the different substrate temperatures. The black line corresponds to the pristine substrates.

Figure

  Figure frequency extracted from the reflectivity spectra as a function of the growth temperature for different substrates.

Figure 8 .

 8 Figure 8. a) c-axis lattice parameter of SVO films grown on STO before (black squares) and after annealing (red stars). b) Absolute variation of c-axis after annealing.

Figure 9 .

 9 Figure 9. Variation of a) the A and b) A 1/2 N -1/3 proportional to the electronic effective mass as a function of growth temperature before and after annealing on STO.

Figure 10 .

 10 Figure 10. Reflectivity spectra of SVO films grown on STO at 400 and 700 °C before (dashed line) and after (solid line) soft annealing at 100 °C-20 h under atmosphere.

44 45 46 47 48 49 50 44 45 46 47 48 49 50 SVO (002) 2 Theta /°3 00°C 400°C 500°C Intensity /a. u. 600°C 700°C STO (002) 44 45 46 47 48 49 50 44 45 46 47 48 49 50 SVO (002) LAO (002) 300°C 2 Theta /°4 00°C 500°C 600°C 700°C Intensity /a. u.

Table 1 .

 1 Transport properties of SVO//STO grown at 400 and 700 °C before and after annealing at 100 °C during 1 day.

	SVO//STO	ρ 300K [× 10 -4 Ω -1 cm -1 ]	ρ 0 [× 10 -4 Ω -1 cm -1 ]	RRR	N [× 10 22 cm 3 ]	Mobility [cm 2 (V -1 s -1 )]
	400 °C As-deposited	1.87	1.19	1.57	2.48	1.35
	400 °C Annealed	1.78	1.14	1.56	2.55	1.37
	700 °C As-deposited	2.62	2.14	1.22	2.09	1.14
	700 °C Annealed	2.07	1.59	1.30	2.26	1.33

Acknowledgements

The authors thank the Labex EMC3 (Energy Materials and Clean Combustion Center) for its financial support in the framework of the COTRA project. The authors thank also M. Boisserie and C. Frilay for important technical help, on the electrical and optical measurements, respectively.

Conflict of Interest

The authors declare no conflict of interest.


Majda Arif

Jérémy Gauthier

Kevin Sugier

Daniele Iudicone

Olivier Jaillon

Patrick Wincker

Pierre Peterlongo

Mohammed-Amin Madoui

Discovering Millions of Plankton Genomic Markers from the Atlantic Ocean and the Mediterranean Sea

Keywords:

Comparison of the molecular diversity in all plankton populations present in geographically distant water columns may allow for a holistic view of the connectivity, isolation and adaptation of organisms in the marine environment. In this context, a large-scale detection and analysis of genomic variants directly in metagenomic data appeared as a powerful strategy for the identification of genetic structures and genes under natural selection in plankton.

Here, we used DiscoSnp++, a reference-free variant caller, to produce genetic variants from large-scale metagenomic data and assessed its accuracy on the copepod Oithona nana in terms of variant calling, allele frequency estimation and population genomic statistics by comparing it to the state-of-the-art method. DiscoSnp++ produces variants leading to similar conclusions regarding the genetic structure and identification of loci under natural selection. DiscoSnp++ was then applied to 120 metagenomic samples from four size fractions, including prokaryotes, protists and zooplankton sampled from 39 Tara Oceans sampling stations located in the Atlantic Ocean and the Mediterranean Sea to produce a new set of marine genomic markers containing more than 19 million of variants. This new genomic resource can be used by the community to relocate these markers on their plankton genomes or transcriptomes of interest. This resource will be updated with new marine expeditions and the increase of metagenomic

Introduction

The identification of population connectivity, isolation and adaptation is of great interest to understand the current and future ecological responses of plankton communities to environmental variations such as the rise of water temperature and acidity [START_REF] Freer | Predicting ecological responses in a changing ocean: the effects of future climate uncertainty[END_REF][START_REF] Pelejero | Paleo-perspectives on ocean acidification[END_REF], especially in a climate change context [START_REF] Beaugrand | Plankton effect on cod recruitment in the North Sea[END_REF][START_REF] Beaugrand | Reorganization of North Atlantic marine copepod biodiversity and climate[END_REF]. To understand the impact of these changes on living organisms, the study of plankton populations at the molecular level is a valuable option since it allows us not only to characterize genetic structures but also to determine which genes and biological functions are under natural selection [START_REF] Avise | Molecular Markers, Natural History and Evolution[END_REF][START_REF] Peijnenburg | High evolutionary potential of marine zooplankton[END_REF]. Previous studies performed on plankton were based mostly on a few molecular markers, such as ribosomal DNA or mitochondrial genes [START_REF] Blanco-Bercial | DNA barcoding of marine copepods: assessment of analytical approaches to species identification[END_REF][START_REF] Cepeda | Molecular systematic of three species of Oithona (Copepoda, Cyclopoida) from the Atlantic Ocean: comparative analysis using 28S rDNA[END_REF]). An alternative capture-based approach based on RADseq has also been proposed [START_REF] Blanco-Bercial | New view of population genetics of zooplankton: RAD-seq analysis reveals population structure of the North Atlantic planktonic copepod Centropages typicus[END_REF]. These approaches permitted the construction of population genetic structures using only a subset of the whole genomic variability. Furthermore, as the loci under selection represent only a very small fraction of a genome, the lack of resolution of these methods does not allow a comprehensive view of the natural selection occurring on plankton. To be able to capture the entire genomic variability of these organisms, whole genome sequencing of individuals could be the ideal strategy. However, due to the small size of certain major zooplankters and their large genome size [START_REF] Wyngaard | Patterns of genome size in the copepoda[END_REF][START_REF] Wyngaard | The relationship between genome size, development rate, and body size in copepods[END_REF], the current DNA extraction methods applied on a single individual do not permit us to retrieve a sufficient amount of genomic DNA that captures the whole genome complexity and that is needed to build genomic DNA libraries (without random genomic amplification) usable for high-throughput sequencing.

Recently, the use of metagenomic data has been proposed to identify natural selection in prokaryotes [START_REF] Costea | metaSNV: A tool for metagenomic strain level analysis[END_REF][START_REF] Delmont | The global biogeography of amino acid variants within a single SAR11 population is governed by natural selection[END_REF][START_REF] Schloissnig | Genomic variation landscape of the human gut microbiome[END_REF]. A similar approach has also been applied to the widespread marine copepod Oithona [START_REF] Madoui | New insights into global biogeography, population structure and natural selection from the genome of the epipelagic copepod Oithona[END_REF] to establish a population genomic analysis at the whole-genome level. The methods used in these studies were all based on metagenomic reads mapping to reference genomes, followed by several filtering steps based on the nucleic identity cut-off and depth of sequencing coverage prior to the variant calling step. This allowed the detection of polymorphic loci and the estimation of allele frequencies in each sample that were followed by a wide range of analyses to characterize the nucleic variations and to identify selection using population genetic metrics such as F ST [START_REF] Wright | The genetical structure of populations[END_REF], LK [START_REF] Lewontin | Distribution of gene frequency as a test of the theory of the selective neutrality of polymorphisms[END_REF], and FLK [START_REF] Bonhomme | Detecting selection in population trees: the Lewontin and Krakauer test extended[END_REF]. In these previous studies, the arbitrary nucleic identity cut-off was used to decrease the amount of false positive variants that can be generated by the alignment of metagenomic reads provided by a closely related species that can be present in the sample. Although the use of such a filter is justified, reads harbouring more variation (< 97% identity) but belonging to the studied organism are de facto discarded. Moreover, the time and computational resources needed for metagenomic read alignments increase with the number of reference genomes included in the analysis. Finally, methods based on read alignments suffer from bias due to the incompleteness and imperfectness of reference genome sequences unless reference genomes are exhaustively and correctly assembled, which is rarely the case.

To bypass these problems, the use of an alignment-free variant calling method could be a solution. Therefore, in the present study, we used DiscoSnp++ [START_REF] Peterlongo | DiscoSnp++: de novo detection of small variants from raw unassembled read set(s) bioRxiv 209965[END_REF][START_REF] Uricaru | Reference-free detection of isolated SNPs[END_REF], a reference-free variant caller, and compared its performance to the one obtained with bwa/samtools/bcftools (BSB) (Li et al. 2009) first using simulated data and then the Tara Oceans metagenomic data [START_REF] Karsenti | A holistic approach to marine eco-systems biology[END_REF][START_REF] Pesant | Open science resources for the discovery and analysis of Tara Oceans data[END_REF] on the O. nana reference genome as a case study to determine DiscoSnp++ accuracy for variant calling, allele frequency estimation and downstream population genomic analysis. Then, we applied DiscoSnp++ to Tara Oceans metagenomic data from the Atlantic Ocean (AO) and the Mediterranean Sea (MS) to provide a new genomic resource that contains more than 19 million marine genomic variants (MGVs) that can be used as is or directly mapped on plankton genomes and transcriptomes of interest for population genomic analysis using the provided DiscoSnp++ module.

Material and Methods

Metagenomic data and genome reference

To compare DiscoSnp++ to BSB, we used metagenomics reads from the MS collected by the Tara Oceans expedition [START_REF] Alberti | Viral to metazoan marine plankton nucleotide sequences from the Tara Oceans expedition[END_REF]) that correspond to the 20-180 µm fraction size from the surface (≤20 m) water layers of Mediterranean stations TARA_8,10,11,12,24 and 26 (Supplementary Notes S1) [START_REF] Pesant | Open science resources for the discovery and analysis of Tara Oceans data[END_REF]. Only the metagenomic data from the two stations TARA_8, and 11 were used to compare the performances of DiscoSnp++ versus BSB for variant calling and B-allele frequency accuracy. Data from the five stations TARA_10,11,12,24 and 26 were used to compare the two approaches in order to perform population genomic analysis. The O. nana genome was downloaded from NCBI (accession number: GCA 900157175.1). To build the marine genomic variants sets (MGVs), we used Tara Oceans metagenomic reads generated from samples corresponding to four size fractions (0.8-5 µm, 5-20 µm, 20-180 µm and 180-2000 µm) collected from stations located in the AO and MS (Supplementary Notes S1).

The BSB pipeline

The bwa mem [START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF] command was used to align the metagenomics reads on the O. nana genome with a 17 bp seed, and alignments were stored in one sorted BAM file per station.

To avoid spurious read alignments, Dust was applied with default parameters to discard reads with low complexity. The reads with an identity under 97% with the O. nana genome were discarded. For the variant (in this study, we will systematically use the term 'variant' to refer to a single nucleotide polymorphism) calling step, we used the samtools mpileup and bcftools call -m commands (Li et al. 2009) with default parameters. Loci with a maximum of two alleles were kept. Only positions with a vertical coverage between the median coverage ± two standard deviations were kept with a minimum of 4x coverage (Supplementary Notes S2).

DiscoSnp++ method overview

DiscoSnp++ was originally designed for genomic data analysis, however, the core of the programme also applies to cases of metagenomic data. The tool is based on the analysis of the de Bruijn Graph (dBG). In the genome assembly context [START_REF] Pevzner | De novo repeat classification and fragment assembly[END_REF]), a dBG is a graph in which nodes are words of length k (k-mers), and each edge connects two k-mers that share a k-1 overlap. For assembling purposes, the dBG is constructed from k-mers of a read set, and contigs are obtained by finding paths in this graph. In practice, k-mers are counted and those having an unexpected low abundance are removed as they are considered to contain sequencing errors. The dBG is constructed with the remaining k-mers. Basically, in a dBG, a bubble denotes a path in the graph which diverges into two distinct paths before they reunite. Any couple of distinct sequences that exists in the data, starting and finishing with the same k-nucleotides, generates a bubble in the dBG. In particular, small indels and SNPs generate such a topological pattern. The DiscoSnp++ algorithm detects bubbles whose couple of paths is of equal length (generated by substitutions in the data) and bubbles whose couple of paths have a difference of length ≤D (generated by insertion or deletion of size at most D). The detection of bubbles in the dBG can be performed through different methods corresponding to different stringencies: parameter -b 0 or 1, with -b 0 providing high precision, and lower recall and conversely (see [START_REF] Uricaru | Reference-free detection of isolated SNPs[END_REF] for more details). In a second step, raw reads are mapped back on the sequence of these paths. This step provides a way to remove non-coherent sequences [START_REF] Myers | The fragment assembly string graph[END_REF] and to supply read coverage per variant and per input read set, whatever the number of input read set(s). This allows the simultaneous analysis of large metagenomic data sets. When a reference genome is available, sequence variants can be mapped to it. Thus, mapped predicted variants have a genomic position, provided in a VCF file.

DiscoSnp++ was run using the default parameters but avoiding indels (-D 0) for the methods comparison, and additionally using -k 51 to build the MGVs. Using a large k value (here k=51) decreases the method sensitivity, increases the precision and, by simplifying the de Bruijn graph, allows faster computing and lower memory use on very large datasets. Depending on the situation, DiscoSnp++ was run using -b 0 (default) or -b 1, and this parameter is specified and motivated in the text. For the BSB pipeline, biallelic loci were kept and only positions with a vertical coverage between the median coverage ± two standard deviations were kept with a minimum of 4x coverage (Supplementary Notes S2).

Comparison of the variant calling methods on simulated data

We simulated a first population of 20 O. nana genomes having 99% identity with the O. nana reference genome and a second population of 20 genomes of "Oithona2" based on a new reference having 95% identity to the O. nana genome and with a 99% identity within the population. SNPs were simulated to reproduce their natural distribution along the genome (Supplementary Notes S3). We generated 100x of Illumina reads on each population and created 20 read datasets by mixing the two populations in different proportions and each dataset contained a total of 30X of simulated Illumina reads. We applied the two approaches to these simulated datasets using the O. nana genome as a reference for reads mapping and variant calling for BSB and variant relocating for DiscoSnp++. The variants found by the two methods (DiscoSnp++ was run in relaxed mode only) were compared to the simulated ones (methodology presented in Figure 1). Considering O. nana as the organism of interest, the signal-to-noise ratio was calculated, i.e., the ratio of the number of O. nana variants over the number of Oithona2 variants.

Comparison of the variant calling methods on real data

The two methods were compared on their performance to identify intra-species variants present in the O. nana genome. The TARA_8 sample was known to contain an abundant species closely related to O. nana with a median identity percentage of 95% and very few O. nana (<10% of total Oithona based on the 28S relative abundance) [START_REF] Madoui | New insights into global biogeography, population structure and natural selection from the genome of the epipelagic copepod Oithona[END_REF]. The variants predicted from this sample by any method and remapped on the O. nana genome can be considered as enriched in inter-species variants. The TARA_11 sample was known to contain a large majority of O. nana (>60% of total Oithona) with other Oithona species that are not closely related to O. nana (here O. similis and O. atlantica). The variants predicted from this sample by any method and remapped on the O. nana genome can be considered enriched in intra-species variants. The stations TARA_8 and 11 were used to compare the two approaches (methodology presented Figure 1), in terms of variant calling, allele frequency accuracy and population genomics statistics (DiscoSnp++ was run in relaxed and stringent mode). To evaluate the possible biases on the coverage of biallelic loci that could be introduced by the variant calling methods, the read depth of the biallelic loci was fitted to a negative binomial distribution and the expected skewness of the distribution was calculated and compared to the observed one. The significance of the method's impact on the coverage skewness was tested by Wilcoxon signedrank tests.

Population genomics analysis

The B-allele frequencies (BAFs), also named alternative allele frequencies compared to a haploid reference genome, were calculated from the VCF files generated by the two methods (DiscoSnp++ in relaxed and stringent mode). Only loci with at least a BAF ≥ 0.05 in one population were selected. To identify populations having the same genomic variant pattern, a PCA was performed based on the BAF of the five populations. To measure the genetic differentiation between the populations, we used the F ST (Wright's fixation index):

𝐹 !" = !(!) !(!)(!!! ! )
, with 𝑝 being a set of BAFs observed in n populations at the same biallelic locus, 𝐸 the mean and 𝑉 the variance. For each locus, the F ST was calculated between each population (pairwise F ST ). The median pairwise F ST was then used to estimate the genetic differentiation between each population.

To evaluate the use of DiscoSnp++ to identify population differentiation, we calculated the median pairwise F ST using four sets of BAFs; set1: BAFs inferred from BSB variants found in common with DiscoSnp++; set2: BAFs inferred from DiscoSnp++ variants found in common with BSB; set3: BAFs inferred from all BSB variants; set4: BAFs inferred from all DiscoSnp++ variants. For the four sets, we used the DiscoSnp++ variants called using the -b 0.

Detection of loci under selection

To detect loci under selection, we calculated the Lewontin-Krakauer (LK) statistic, which is an improvement of the F ST that can be used for testing the neutrality of polymorphic genes,

𝐿𝐾 = (!!!)! !"
!(! !" ) . To be able to detect loci under selection, the LK distribution must follow a chisquare distribution χ 2 (n -1) with n being the number of different populations. The fitting between the theoretical χ 2 distribution and the observed LK distribution obtained from BSB and DiscoSnp++ (with -b 0 and -b 1 options) was observed to validate the neutral model, i.e., the majority of the biallelic loci are not under selection (Supplementary Notes S4). The FLK statistics [START_REF] Bonhomme | Detecting selection in population trees: the Lewontin and Krakauer test extended[END_REF] were also calculated; this metric is an extension of the LK test that uses a kinship matrix of the populations based on the BAF to correct genetic distance biases due to population structure. The FLK statistics were also tested for the neutral model. The first hundred loci having the highest LK or FLK values higher than expected (with a p-value ≤ 0.05) were considered to be under selection. The annotation of the variants and their possible effect on protein structure was performed with SnpEff [START_REF] Cingolani | A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3[END_REF].

Results

Variant calling

The BSB and DiscoSnp++ pipelines were compared for variant detection (methodology presented in Figure 1) using simulated data representing an admixture of O. nana and a closely related species in different proportions. Here, we considered O. nana as the organism of interest, its variants were considered as true positives and the variants of Oithona2 as false positives. BSB found more true positives than DiscoSnp++ in all admixtures especially for low O. nana content (between 5 and 50% of O. nana) (Figure 2.a, Supplementary Notes S5). BSB also identified more false positives than DiscoSnp++ especially for admixtures with O. nana lower than 90%. Based on these simulations, DiscoSnp++ was less sensitive in any admixture but more specific than BSB when dealing with an admixture of two closely related species. However, the signal-to-noise ratio was higher for BSB, especially for admixtures with more than 90% of O. nana (Figure 2.b, Supplementary Notes S5).

Allele frequency accuracy

The BAFs obtained for variants found by the two calling methods were compared (Figure 4) and we observed a strong correlation between the two methods in O. nana populations from TARA_8 and 11 (R 2 ≥ 0.95). However, we found that 7.5% of the variants had a higher BAF difference than expected between the two methods (i.e., with a BAF difference higher/lower than the median difference plus/minus two standard deviations, see Supplementary Notes S6.a). For variants having a higher BAF with DiscoSnp++ (6.3% of the total variants found in common), we explained the difference by the identity cut-off of 97% used in the BSB pipeline (Supplementary Notes S7). The variants presenting a strong BAF deviation between the two methods were annotated based on their genomic location (i.e., intronic, exonic, UTR and intergenic) and compared to (i) the genomic location of the variants presenting no significant BAF differences, and (ii) a random distribution of the variants on the genome. Significant differences (p-value < 0.001, chi-square test) were found, with an increase of biallelic loci having higher BAFs with DiscoSnp++ located in the non-coding regions of the genome (Supplementary Notes S6.b and c). This result suggests that DiscoSnp++ can recruit more reads than BSB in noncoding regions of the genome. These regions are indeed expected to contain more polymorphisms than coding regions within populations. Therefore, filters that are applied in BSB tend to discard reads that should be aligned at a reduced similarity threshold. Consequently, DiscoSnp++ seems to provide a better estimation of the allele frequency in more variable regions of the genome compared to BSB applied with a 97% identity cut-off.

Population genomic analysis

Five O. nana populations from sampling stations (TARA_10, 11, 12, 24 and 26) were clustered by PCA based on their BAFs (Figure 5.A and 5.B). For the two methods, the clustering showed a 1).

Using only variants detected by the two methods (i.e., using BAF from set 1 and set 2), we found no F ST difference over 0.024 and the average difference between the median pairwise F ST was 0.012 ± 0.01 (Figure 5.C). A higher difference was observed for high F ST values. We found a negligible difference between pairwise F ST computed by DiscoSnp++ in relaxed mode versus stringent mode (mean=0.003, sd=0.008). For the selected variants, the two methods allowed the identification of the same genetic pattern between the five O. nana populations. The genetic distance observed using all DiscoSnp++ or all BSB variants also produced a similar genetic pattern with an absence of genetic structure within the WMS and a weak differentiation between the two MS basins and within the EMS (Table 1). Compared to previously published results [START_REF] Madoui | New insights into global biogeography, population structure and natural selection from the genome of the epipelagic copepod Oithona[END_REF], there was a lower genetic distance between the population of TARA_26 and the four other populations. This difference can be due to the more stringent filtering on reads coverage used in the current study (see Materials and Methods) to consider valid variants compared to the previous study where biallelic loci with a read coverage up to 80x were kept.

The current coverage filters may have discarded reads provided by repeated regions or a closely related species possibly present in the TARA_26 sample.

Detection of loci under natural selection

To identify loci under natural selection, the LK and FLK statistics were computed from the BAFs of sets 1 and 2. For each variant set and statistics, the hundred loci with the highest LK and FLK were compared to estimate the congruence between the two variant calling methods (Figure 6).

We found more loci in common with LK than FLK and by using the -b 0 option of DiscoSnp++ suggesting a more accurate detection of loci under selection as being more stringent in the variant calling of DiscoSnp++.

The functional annotation of the 79 variants under natural selection (detected by DiscoSnp++ -b 0 and using the LK outliers) that were found in common with BSB (Supplementary Notes S7)

showed 16 non-synonymous variants and 14 synonymous variants. Compared to the previous study [START_REF] Madoui | New insights into global biogeography, population structure and natural selection from the genome of the epipelagic copepod Oithona[END_REF], we found four new Lin12 Notch Repeat (LNR) domain-coding genes.

These domain-coding genes are of particular interest in O. nana where they were found to be over-abundant compared to other metazoans and one of them detected under positive selection was male-specific based on expression data [START_REF] Madoui | New insights into global biogeography, population structure and natural selection from the genome of the epipelagic copepod Oithona[END_REF]). Among the four new LNR

Plankton genomic variant resources from the Tara Oceans metagenomic data

We produced the new set of MGVs by running DiscoSnp++ in relaxed mode (to optimize the number of MGVs) on more than 40 billion metagenomic 100 bp reads from 39 Tara stations located in the AO and the MS (Figure 7 and Supplementary Notes S8). These MGVs correspond to genomic variants (SNVs and indels) found from natural populations of prokaryotic, protist and animal plankton that were sampled during the three-year expedition of Tara. For the four different size fractions, we generated more than nineteen million MGVs (Table 2). The amount of input data was relatively similar among all size fractions (~11-12.10 9 of 100 bp reads) but the computation time globally increased with the size fraction and all had the same very low memory usage (~100 Gb). The amount of MGVs found in the different fraction sizes was at the same scale (5.2-6.2.10 6 variants) except for fraction 5-20 µm that presented half the MGVs and had the lowest computation time. This may be because of less genomic complexity in this fraction size, as shown previously [START_REF] Carradec | A global ocean atlas of eukaryotic genes[END_REF]. The MGVs can be downloaded and directly used by the scientific community in order to perform new analyses of genomic diversity on any organism of interest as demonstrated on O. nana.

Discussion

Like any reference based variant detection method, DiscoSnp++ limitations are mainly due to genomic approximate repeats. Reads from approximate repeats and, in the metagenomic framework, reads from similar inter-species genomic regions contain the same signal as those from regions containing intra-species variants. As shown by the results from this study, those imperfect predictions are not an insurmountable limitation for population genomics analysis where alignment-based and reference-free-based approaches provide similar conclusions in terms of population differentiation and overlapping results in terms of natural selection. Moreover,

Figure 1 :

 1 Figure 1: Workflow for BSB and DiscoSnp++ methods comparison.

Figure 2 :Figure 3 :

 23 Figure 2: Comparison of variant calling between DiscoSnp++ and BSB on simulated data. a. Variants recall for increasing proportion of O. nana in the admixture. b. Methods efficiency for increasing proportion of O. nana in the admixture.

Figure 4 :

 4 Figure 4: B-Allele frequency correlation between DiscoSnp++ and BSB. The x-axis of the scatter-plots corresponds to BAFs obtained with BSB and the y-axis corresponds to BAFs obtained with DiscoSnp++. a. y-axis is BAFs computed from TARA_8 with DiscoSnp++ -b 0. b. y-axis is BAFs computed from TARA_8 with DiscoSnp++ option -b 1. c. y-axis is BAFs computed from TARA_11 with DiscoSnp++ -b 0. d. y-axis is BAFs computed from TARA_11 with DiscoSnp++ -b 1.

Figure 5 :

 5 Figure 5: O. nana genetic structure in the Mediterranean Sea obtained with DiscoSnp++ and BSB. a. PCA on five O. nana populations on the MS based on BAFs obtained with BSB. b. PCA on five O. nana populations on the MS based on BAFs obtained with DiscoSnp++. c. Differences of the median pairwise F ST between BSB and DiscoSnp++.

Figure 6 :

 6 Figure 6: Loci under natural selection found in common between DiscoSnp++ and BSB. a. Loci in common using the 100 LK highest values (p-value<0.001, chi-square test). b. Loci in common using the 100 highest FLK values (p-value<0.001, chi-square test).

  domain-coding genes found to be under selection, one (GSONAT00015400001) codes a metallopeptidase domain protein, another (GSONAT00015380001) codes an LNR protein associated with a Kelch domain and two others (GSONAT00013822001, GSONAT00015410001) code only LNR domain proteins without association to other known domains. These new results reinforce the highly evolutionary potential of LNR domaincontaining proteins and their importance in the O. nana biology.

Figure 7 :

 7 Figure 7: Geographic and size fraction distribution of MGVs.

Table 1 : Median pairwise F ST between O. nana populations obtained from the four BAFs sets.

 1 Set1: BAFs inferred from BSB for all variants found in common with DiscoSnp++; set2: BAFs inferred from DiscoSnp++ for variants found in common with BSB; set3: all BAFs inferred from BSB; set4: all BAFs inferred from DiscoSnp++. grouping of the populations by geographic location, separating the ones from the Western MS (WMS) from the ones of the Eastern MS (EMS). We estimated the genetic differentiation between the O. nana populations by calculating the pairwise F ST using four different sets of variants and related BAFs (Supplementary Notes S6) and compared the median F ST values to evaluate any biases that could be introduced by DiscoSnp++ ran in stringent mode (Table

			Median pairwise F ST
						Standard
	Populations	Set 1	Set 2	Set 3	Set 4	deviation
	TARA_10 vs TARA_11	0.074	0.074	0.065	0.075	0.0046
	TARA_10 vs TARA_12	0.077	0.084	0.065	0.086	0.0096
	TARA_10 vs TARA_24	0.096	0.096	0.077	0.099	0.01
	TARA_10 vs TARA_26	0.109	0.125	0.096	0.133	0.016
	TARA_11 vs TARA_12	0.077	0.099	0.071	0.1	0.0149
	TARA_11 vs TARA_24	0.100	0.124	0.089	0.128	0.0189
	TARA_11 vs TARA_26	0.114	0.142	0.099	0.143	0.0216
	TARA_12 vs TARA_24	0.096	0.1	0.077	0.105	0.0121
	TARA_12 vs TARA_26	0.096	0.1	0.077	0.111	0.0141
	TARA_24 vs TARA_26	0.105	0.116	0.095	0.125	0.013

similar

Table 2 : Marine genomic variants produced by DiscoSnp++ on Tara Oceans metagenomic data from the Atlantic Ocean and the Mediterranean Sea.

 2

	Fraction size	Number of	Number of	Number of	Computation	Max memory
	(μm)	stations	reads used	variants	time (hours)	used (Gb)
	0.8-5	25	11.3 x 10 9	5.5 x 10 6	64	107
	5-20	27	11.8 x 10 9	2.3 x 10 6	60	107
	20-180	31	11.2 x 10 9	5.2 x 10 6	105	110
	180-2 000	31	11.2 x 10 9	6.2 x 10 6	124	120

Acknowledgements

We wish to thank the individuals and sponsors who participated in the Tara Oceans Expedition 2009-2013: Centre National de la Recherche Scientifique, European Molecular Biology Laboratory, Genoscope/Commissariat à l'Energie Atomique, the French Government "Investissements d'Avenir" programmes OCEANOMICS (ANR-11-BTBR-0008), FRANCE GENOMIQUE (ANR-10-INBS-09-08) and HYDROGEN (ANR-14-CE23-0001). This is Tara Oceans contribution number 83.

BWA/samtools/bcftools (BSB), Mediterranean Sea (MS), Atlantic Ocean (AO)

http://bioinformatique.rennes.inria.fr/taravariants/

Data availability

The metagenomic data from Tara Oceans are available at ENA (Supplementary data S1). The Oithona nana genome sequence and annotation are available at ENA with the study Accession no. PRJEB18938. The MGVs files and their corresponding tutorial are available at http://bioinformatique.rennes.inria.fr/taravariants/.

DiscoSnp++ is an order of magnitude faster and uses fewer resources [START_REF] Peterlongo | DiscoSnp++: de novo detection of small variants from raw unassembled read set(s) bioRxiv 209965[END_REF][START_REF] Uricaru | Reference-free detection of isolated SNPs[END_REF]. In the case of the admixture of two closely related species, neither the alignment-based nor the reference free-based approach allows the removal of inter-species variants which reduce the number of populations that can be integrated into the population genomic analysis focused on a single species.

The MGVs detected de novo with DiscoSnp++ from Tara Oceans data can be downloaded from http://bioinformatique.rennes.inria.fr/taravariants/ and used directly on any genome or transcriptome provided by the users to create VCF files without computation of the variant calling. This can be done by running only the final DiscoSnp++ step 'run_VCF_creator.sh' that can be done on a laptop computer. This allows the community to avoid (i) the systematic downloading of the whole Tara Oceans metagenomic data set that needs investment in large infrastructure for data storage and backup and (ii) the alignment of the reads to their genomes and transcriptomes of interest that needs investment in computational power. As demonstrated in this study, the MGVs allow an accurate analysis of the molecular diversity of the plankton present in the AO and MS that were captured during the Tara Oceans expedition. In addition to the lack of reference sequences for plankton, depending on the genome size and abundance of the studied plankton in the Tara Oceans samples, the use of the MGVs collection may have some limits.

Analyses focusing on small-size genomes (<100 Mb) and abundant protists such as green algae are more likely to provide interesting results compared to those focusing on copepods with largesize genomes (>1 Gb).

The increasing number of large collections of marine plankton samples and their related metagenomic dataset forces a rethinking of the way population genomics can be performed. This can push the community towards the use of a universal genomic resource of variants that can be updated with the accumulation of newly released metagenomic data. From this perspective, the use of DiscoSnp++ offers a great advantage by providing a uniform method to generate community shared markers that store all the information needed to perform robust downstream population genetics analyses of plankton.

Authors' contribution

MA, KS, PP, JG and MAM performed the analyses. OJ, DL, PP and MAM designed the study.

MAM wrote the manuscript, and all authors accepted its final version.