Jintai Fan

Vincent Sarou-Kanian

Xiaoyan Yang

Maria Diaz-Lopez

Franck Fayon

Xiaojun Kuang

Michael J Pitcher

Mathieu Allix

a metastable ternary melilite with a super

ABSTRACT:

The La1+xAE1-xGa3O7+x/2 melilite family (AE = Ca, Sr, Ba, and 0  x  0.64) demonstrates remarkable oxide ion conductivity due to the ability of its layered tetrahedral [Ga3O7+x/2] network to accommodate and transport interstitial oxide ions (Oint). Compositions x > 0.65 with very high Oint concentrations (referred to here as "super-excess" compositions) have the potential to support correspondingly high ionic conductivities, but have never before been accessed due to the limitations of conventional solidstate ceramic synthesis. Here we report that fully-substituted La2Ga3O7.5 (x = 1) melilite ceramics can be synthesized by direct crystallization of an under-cooled melt, demonstrating that super-excess compositions are accessible under suitable nonequilibrium reaction conditions. La2Ga3O7.5 is stable up to 830 °C and exhibits an ionic conductivity of 0.01 S.cm -1 at 800 °C, three orders of magnitude higher than the corresponding x = 0 end-member LaSrGa3O7, and close to the range exhibited by the current best-in-class La1.54Sr0.46Ga3O7.23 (0.1 S.cm -1 ). It crystallizes in an orthorhombic √2a x √2a x 2c expansion of the parent melilite cell, in space group Ima2, with full long-range ordering of Oint into chains within the [Ga3O7.5] layers. The emergence of this chain-like (1D) ordering within the 2D melilite framework, which appears to be an incipient feature of previously reported partially-ordered melilites, is explained in terms of the underlying hexagonal topology of the structure. These results will enable the exploration of extended compositional ranges for the development of new solid oxide ion electrolytes with high concentrations of interstitial oxide charge carriers.

The discovery of new compounds with high oxide ion conductivities is important for the development of solid oxide fuel cell (SOFC) technology, where it is desirable to reduce device operating temperatures to the range 500 -800°C to improve durability and reduce running costs. [START_REF] Wachsman | Lowering the Temperature of Solid Oxide Fuel Cells[END_REF][START_REF] Kilner | Materials for Intermediate-Temperature Solid-Oxide Fuel Cells[END_REF] However, it can be challenging to access new compositions using classic ceramic synthesis routes, because the high temperatures and extended heating times required tend to favour the formation of thermodynamic products. This means that potential metastable compositions and structures, which may be highly stable at the temperatures of interest for their function, can be missed. Aerodynamic levitation (ADL) coupled to laser heating [START_REF] Winborne | Aerodynamic Levitation Technique for Containerless High Temperature Studies on Liquid and Solid Samples[END_REF][START_REF] Benmore | Aerodynamic Levitation, Supercooled Liquids and Glass Formation[END_REF] offers routes to new metastable bulk oxides by the crystallization-from-glass method, where a glass precursor is fabricated by rapid quenching of a congruent melt and then crystallised by annealing above the glass transition temperature, or by direct crystallisation of the under-cooled melt. Both routes have recently emerged as useful processes for the non-equilibrium synthesis of new bulk oxides in the form of high-density low-porosity ceramics, [START_REF] Allix | Highly Transparent BaAl4O7 Polycrystalline Ceramic Obtained by Full Crystallization from Glass[END_REF][START_REF] Véron | Synthesis and Structure Determination of CaSi1/3B2/3O8/3: A New Calcium Borosilicate[END_REF][START_REF] Boyer | BaGa4O7, a New A3BC10O20 Crystalline Phase: Synthesis, Structural Determination and Luminescence Properties[END_REF][START_REF] Boyer | First Oxide Ion Conducting Ceramics Synthesized by Full Crystallization from Glass[END_REF][START_REF] Ma | Pressureless Glass Crystallization of Transparent Yttrium Aluminum Garnet-Based Nanoceramics[END_REF] and are amenable to chemistries that support interstitial oxide ion conductivity, as demonstrated recently for gallate- [START_REF] Boyer | First Oxide Ion Conducting Ceramics Synthesized by Full Crystallization from Glass[END_REF] and aluminate-based [START_REF] Xu | Interstitial Oxide Ion Migration Mechanism in Aluminate Melilite La1+xCa1-xAl3O7+0.5x Ceramics Synthesized by Glass Crystallization[END_REF] melilites and the scheelite-type bismuth vanadates. [START_REF] Yang | Cooperative Mechanisms of Oxygen Vacancy Stabilization and Migration in the Isolated Tetrahedral Anion Scheelite Structure[END_REF] The melilite structure has been of interest for potential SOFC applications since the first reports of high ionic conductivity in non-stoichiometric La1+xSr1-xGa3O7+x/2. [START_REF] Rozumek | Electrical Conduction Behavior of La1+xSr1-xGa3O7-δ Melilite-Type Ceramics[END_REF] Its structure (general formula A2B3O7) consists of a two-dimensional framework of corner-sharing 3-and 4-connected BO4 tetrahedra forming distorted pentagonal channels perpendicular to the stacking axis, which alternates with a layer of A cations in 8-fold coordination by oxide ions. In La1+xAE1-xGa3O7+x/2 (AE = Ca 2+ , Sr 2+ or Ba 2+ ), interstitial oxide ions (Oint) are incorporated into the framework close to the centres of the pentagonal rings by tuning the La 3+ /AE 2+ ratio (x). [START_REF] Zhou | Development of Melilite-Type Oxide Ion Conductors[END_REF] This creates GaO5 trigonal bipyramidal units locally by coordination of Oint to gallium ions within the tetrahedral layer, which is stabilized by a substantial bonding contribution from the two adjacent La 3+ /Sr 2+ cations. [START_REF] Kuang | Interstitial Oxide Ion Conductivity in the Layered Tetrahedral Network Melilite Structure[END_REF][START_REF] Xu | Chemical Bonding Effect on the Incorporation and Conduction of Interstitial Oxide Ions in Gallate Melilites[END_REF] At elevated temperatures, these interstitial oxide ions act as chargecarriers, primarily by diffusion within the Ga-O framework layer. [START_REF] Tealdi | Energy Materials: Layered LaSrGa3O7-Based Oxide-Ion Conductors: Cooperative Transport Mechanisms and Flexible Structures[END_REF][START_REF] Wei | Anisotropic Oxide Ion Conduction in Melilite Intermediate Temperature Electrolytes[END_REF] Consequently, the ionic conductivity (σ) depends directly on the concentration of interstitial oxide (x) in the framework. In La1+xAE1-xGa3O7+x/2 (AE = Ca 2+ , Sr 2+ ), σ shows a strong increase with x, and reaches a broad maximum of ~0.1 Scm -1 at 800°C when 0.5 ≤ x ≤ 0.6. [START_REF] Wei | Anisotropic Oxide Ion Conduction in Melilite Intermediate Temperature Electrolytes[END_REF][START_REF] Li | Interstitial Oxide Ion Order and Conductivity in La1.64Ca0.36Ga3O7.32 Melilite[END_REF][START_REF] Thomas | Phase Stability Control of Interstitial Oxide Ion Conductivity in the La1+xSr1-xGa3O7+x/2 Melilite Family[END_REF] At higher concentrations (0.60 < x ≤ 0.65 for AE = Sr 2+ and 0.55 ≤ x ≤ 0.64 for AE = Ca 2+ ), the onset of long-range oxide ordering induces an orthorhombic distortion, which is associated with a slight decrease in σ. [START_REF] Thomas | Phase Stability Control of Interstitial Oxide Ion Conductivity in the La1+xSr1-xGa3O7+x/2 Melilite Family[END_REF] The compositional limit of the solid solution is finally reached at x = 0.65: above this limit, single-phase melilite samples cannot be obtained from AE = Ca 2+ , Sr 2+ by ceramic synthesis. [START_REF] Li | Interstitial Oxide Ion Order and Conductivity in La1.64Ca0.36Ga3O7.32 Melilite[END_REF][START_REF] Thomas | Phase Stability Control of Interstitial Oxide Ion Conductivity in the La1+xSr1-xGa3O7+x/2 Melilite Family[END_REF] A similar trend in ionic conductivity with increasing x is observed in La1+xBa1-xGa3O7+x/2, which has a smaller compositional stability range of 0 ≤ x < 0.35 and remains tetragonal throughout the series. [START_REF] Xu | La1+xBa1-xGa3O7+0.5x Oxide Ion Conductor: Cationic Size Effect on the Interstitial Oxide Ion Conductivity in Gallate Melilites[END_REF] In all three series, the maximum observed interstitial oxide concentration (x = 0.65) is substantially lower than the proposed saturation value of the framework at x = 1, based on the consideration of local strain effects in the Ga-O layer, [START_REF] Li | Interstitial Oxide Ion Order and Conductivity in La1.64Ca0.36Ga3O7.32 Melilite[END_REF] which would correspond to full substitution of AE by La 3+ . This implies that the structure may be amenable to further enrichment with interstitial oxide if suitable thermodynamic routes can be identified. The ability to break through this limit to super-excess Oint concentrations 0.65 < x ≤ 1 would substantially increase the explorable phase space in this class of compounds, raising the prospect of further improvements in conductivity and new structural chemistry. Recent studies suggest that non-equilibrium synthesis methods can be used to exceed such compositional limits. For example, a sol-gel approach with annealing at moderate temperatures drastically extends the accessible interstitial oxide concentration range in langasite-type La3Ga5-xGe1+xO14+x/2; 21 whilst the synthesis of metastable melilites including the nonstoichiometric aluminates La1+xCa1-xAl3O7+x/2 (x ≤ 0.6) and a range of stoichiometric and nonstoichiometric small-A-site compositions, which cannot be achieved by high temperature ceramic routes, is made possible by crystallization-fromglass. [START_REF] Boyer | First Oxide Ion Conducting Ceramics Synthesized by Full Crystallization from Glass[END_REF][START_REF] Xu | Interstitial Oxide Ion Migration Mechanism in Aluminate Melilite La1+xCa1-xAl3O7+0.5x Ceramics Synthesized by Glass Crystallization[END_REF][START_REF] Boyer | Transparent Polycrystalline SrREGa3O7 Melilite Ceramics: Potential Phosphors for Tuneable Solid State Lighting[END_REF] A recent study of (1-y)La2O3 -(y)Ga2O3 glasses by ADL-laser heating found that y = 0.40 was difficult to vitrify despite lying well within the limits of the glass-forming region (0.20 ≤ y ≤ 0.55), hinting at the presence of an unknown metastable compound. [START_REF] Yoshimoto | Low Phonon Energies and Wideband Optical Windows of La2O3-Ga2O3 Glasses Prepared Using an Aerodynamic Levitation Technique[END_REF] Motivated by the potential for a non-equilibrium route to super-excess melilites, we describe here the synthesis of the previously unknown ternary compound La2Ga3O7.5 (y = 0.40) and its structural and physical characterisation. We found that this compound forms in a single step by rapid direct crystallisation of the under-cooled melt, and adopts a new fully-ordered melilite superstructure with a chain-like arrangement of interstitial oxide ions. La2Ga3O7.5, the x = 1 end member of the La1+xAE1-xGa3O7+x/2 interstitial oxide conductor family, is stable up to 830°C and retains a high oxide ion conductivity, suggesting that long-range oxide ion ordering is only moderately detrimental to electrolyte performance and opening the door to the exploration of a new range of melilite compositions. EXPERIMENTAL Synthesis procedure. High-purity La2O3 powder (STREM 99.99%, pre-dried at 1100 °C for 12 hours in air) and Ga2O3 powder (STREM, 99.999%) were weighed according to the molar ratio 39La2O3-61Ga2O3 (where the excess of Ga2O3 compared to the target composition La2Ga3O7.5 was used to compensate for vaporization during the laser heating steps). The La2O3 and Ga2O3 powders were mixed in an agate mortar under ethanol, then dried at 100 °C to remove residual ethanol. The dry powder was then pressed into pellets under 20 MPa uniaxial pressure. Pellet fragments of ~30 mg were placed into the conical nozzle of the aerodynamic levitator system. [START_REF] Winborne | Aerodynamic Levitation Technique for Containerless High Temperature Studies on Liquid and Solid Samples[END_REF][START_REF] Weber | The Containerless Synthesis of Glass[END_REF][START_REF] Allix | Aerodynamic levitation[END_REF] The samples reported in this work were heated by a pair of CO2 lasers while levitating in a flow of oxygen gas, with temperature mon-itored in-situ by an IR pyrometer. Melting of the reaction mixture occurred at approximately 1650 °C, and the samples were further heated up to 1800 °C for several seconds to homogenize the melt. The temperature was then decreased to 1500 °C before switching off the lasers, thus quenching the samples to room temperature. The occurrence of direct crystallization of the melt was indicated by an intense flash from the cooling bead. This heating and quenching process was repeated several times until the beads adopted a homogeneous transparent appearance (see Figure S1). No further heat treatment was applied to the quenched beads. The beads were then either crushed and ground into powders for structural characterization, or were optically polished with an automatic polisher on silicon carbide papers to a thickness of approximately 1 mm for conductivity measurements. Our attempts to synthesize high quality samples of La2Ga3O7.5 by low temperature glass-crystallization were unsuccessful due to the difficulty in vitrifying compositions close to 39La2O3-61Ga2O3. Further details, including the crystallization behavior of the strongly Ga-deficient glasses that we obtained (Figures S2 andS3), are provided in the Supporting Information. Powder Diffraction. Provisional phase identification and indexing was conducted using a Bruker D8 Advance diffractometer with Cu Kα radiation in Bragg-Brentano mode. Variable temperature PXRD (VT-PXRD) data were recorded during heating from room temperature up to 850 °C under air using a HTK1200N Anton Paar furnace. The powder sample was placed in a platinum-lined corundum sample holder. Data were collected from 12 -70° (2θ) with a 0.016 ° step size and an acquisition time of 0.5s per step. The evolution of the unit cell parameters was obtained from Le Bail refinements. For combined synchrotron powder X-ray diffraction (SPD) and neutron powder diffraction (NPD) experiments a large (~5 g) powder sample of La2Ga3O7.5 was prepared by crushing and combining a large number of individual beads. SPD data were collected at 11-BM at the Advanced Photon Source, Argonne National Laboratory, U.S.A. The sample was contained in a Kapton capillary of inner diameter 0.8mm. Data were collected from a spinning sample (60 Hz) over the 0.5-50° 2θ range with a 0.001° step size and acquisition duration of 0.1 s per step at room temperature with λ = 0.412737 Å. NPD experiments were carried out at the D2B high-resolution/high-flux powder diffractometer at the Institut Laüe-Langevin in Grenoble, France. The samples were packed in vanadium cans and data were acquired with λ = 1.594 Å at 2θ intervals of 0.05° over the 10°< 2θ < 160° angular range. The structure of the sample was refined using a combination of SPD and NPD data using the Rietveld method implemented in Jana2006. [START_REF] Petrícek | Crystallographic Computing System JANA2006: General Features[END_REF] Bond valence sum maps were calculated using BondStr. [START_REF] Rodrigues-Carvajal | FullProf Suite[END_REF] X-ray total scattering data were collected at the I15-1 beamline at the Diamond Light Source, UK, using λ = 0.161669 Å (76.7 keV). A sample of La2Ga3O7.5 was contained in a borosilicate glass capillary with an inner diameter of 1 mm and measured for several minutes at room temperature using a Perkin Elmer area detector placed ~200mm downstream from the sample. The position of the sample was calibrated with a Si standard and the collected 2D X-ray scattering data was processed into 1D patterns using the DAWN software. [START_REF] Filik | Processing Two-Dimensional X-Ray Diffraction and Small-Angle Scattering Data in DAWN 2[END_REF] The data were corrected using equivalent measurements taken from an empty capillary and Fourier-Transformed using the GudrunX software with Qmax = 25Å -1 . A simultaneous Rietveld and pair distribution function (PDF) refinement using the same model was performed with TOPAS v.6. [START_REF] Coelho | TOPAS Academic version 6[END_REF] An r-dependent function of the thermal parameters as implemented in TOPAS v.6 was used to fit the data within the range 1.5-50 Å, where the Biso of all atoms were constrained to match those in Rietveld for r ≥ 20 Å. Calculations. Density Functional Theory (DFT) calculations with periodic boundary conditions were performed with the CASTEP code [START_REF] Segall | First-Principles Simulation: Ideas, Illustrations and the CASTEP Code[END_REF][START_REF] Clark | First Principles Methods Using CASTEP[END_REF] using the Materials Studio 7.0 interface. The PBE functional [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] was used to describe electron exchange-correlation effects. Core-valence interactions were described by ultrasoft pseudopotentials generated on-the-fly. [START_REF] Yates | Calculation of NMR Chemical Shifts for Extended Systems Using Ultrasoft Pseudopotentials[END_REF] Geometry optimization was performed using the BFGS minimization method [START_REF] Pfrommer | Relaxation of Crystals with the Quasi-Newton Method[END_REF] and fixing the cell parameters to experimental values. Electric field gradient (EFG) and NMR chemical shift tensors computations were carried out for the experimentally refined structural models using the projector-augmented waves (PAW) [START_REF] Profeta | Accurate First Principles Prediction of 17O NMR Parameters in SiO2: Assignment of the Zeolite Ferrierite Spectrum[END_REF] and gauge included PAW (GIPAW) [START_REF] Pickard | All-Electron Magnetic Response with Pseudopotentials: NMR Chemical Shifts[END_REF] algorithms, respectively. A kinetic energy cutoff of 600 eV was used for the plane wave basis set expansion and the Brillouin zone was sampled using a Monkhorst-Pack grid spacing of 0.04 Å -1 (corresponding to a k-point mesh of 3×3×3). The 71 Ga isotropic chemical shifts (δISO) were deduced from calculated isotropic shielding values (σISO) using the relation δISO = -0.912 σISO + 1519, which accounts for systematic errors on computed values using this method. This relationship was obtained from the linear regression between experimental δISO and calculated σISO values for seven crystalline gallium-containing compounds of known structures (see Table S1 and Figure S5). For all reference compounds, computations were performed without DFT-geometry optimization using the same parameters as those described above. For structural models, 1×2×1 supercells were used (containing 200 atoms). DFT optimization of the atomic positions and cell parameters and total energy computations were performed using ultrasoft pseudopotentials requiring lower energy cutoff of 380 eV with a Monkhorst-Pack grid spacing of 0.04 Å -1 (corresponding to a k-point mesh of 2x1x2). Solid-state Nuclear Magnetic Resonance. We selected 71 Ga (nuclear spin I = 3/2) as the NMR-active nucleus due to the large quadrupolar broadenings of 69 Ga (I = 3/2) and 139 La (I = 7/2), and the low abundance of [START_REF] Wei | Anisotropic Oxide Ion Conduction in Melilite Intermediate Temperature Electrolytes[END_REF] O (I = 5/2). 71 Ga solid-state NMR experiments were performed on a Bruker Avance Neo spectrometer operating at 20.0 T ( 71 Ga Larmor frequency of 259.3 MHz) with a 0.7 mm magic angle sample spinning (MAS) double-resonance probehead. The MAS frequency was set to 100 kHz to avoid overlap between the central transition (CT) centre band and its associated spinning sidebands. The 1D quantitative CT-selective 71 Ga spectrum was recorded using a Hahn-echo experiment with acquisition at the echo top using an inter-pulse delay of 20 µs. The radio frequency (RF) field was 125 kHz to ensure a selective and homogenous excitation of the whole CT spectral range (90-and 180-degree pulse durations of 1.0 and 2.0 µs, respectively). The 2D double-quantum filtered (DQF-) STMAS spectrum [START_REF] Gan | Isotropic NMR Spectra of Half-Integer Quadrupolar Nuclei Using Satellite Transitions and Magic-Angle Spinning[END_REF][START_REF] Kwak | Double-Quantum Filtered STMAS[END_REF] was recorded using an RF field of 420 kHz for the excitation and mixing pulses with pulse durations of 0.4 µs and 0.3 µs, respectively. For the DQ filter and readout pulses the RF-field was 125 kHz. 124 rotor synchronized t1-increments with 4864 transients each were recorded with a recycle delay of 0.25 µs. The magic angle was carefully set by first performing a [START_REF] Yoshimoto | Low Phonon Energies and Wideband Optical Windows of La2O3-Ga2O3 Glasses Prepared Using an Aerodynamic Levitation Technique[END_REF] Na 2D DQF-STMAS experiment on a polycrystalline NaH2PO4 sample. 71 Ga chemical shifts were referenced relative to a 1.1 M aqueous solution of Ga(NO3)2. All spectra were reconstructed using the Dmfit software. [START_REF] Massiot | Modelling One-and Two-Dimensional Solid-State NMR Spectra[END_REF] AC Impedance Spectroscopy. Impedance spectroscopy measurements were carried out with a Solartron 1260 frequency response analyzer over the 10 -1 -10 7 Hz frequency range within 300-800 °C. Platinum paste was coated on the opposite faces of the polished disks and fired at 600 °C for 30 min to burn out the organic components to form electrodes. The impedance data analysis was carried out using the ZView software.

RESULTS AND DISCUSSION

Synthesis and Crystal Structure. Laboratory PXRD on crushed beads of La2Ga3O7.5 synthesised by direct crystallisation of the under-cooled melt produced a new diffraction pattern that did not correspond to any mixture of known compounds, but closely resembled that of the pseudo-orthorhombic melilite La1.64Ca0.36Ga3O7.32 reported by Li et. al. [START_REF] Li | Interstitial Oxide Ion Order and Conductivity in La1.64Ca0.36Ga3O7.32 Melilite[END_REF] (see Figure S1). Compositional analysis of a polished crystalline bead by SEM-EDX showed a Ga/La ratio of 1.51(1), consistent with La2Ga3O7.5. Together, these results suggest that the crystalline beads of La2Ga3O7.5 contain a new phase whose structure is derived from melilite. Structural analysis was performed using SPD data and constantwavelength neutron diffraction (NPD) data simultaneously. To generate a provisional structural model for La2Ga3O7.5, we used the published pseudo-orthorhombic (P1) structure of La1.64Ca0.36Ga3O7.32 [START_REF] Li | Interstitial Oxide Ion Order and Conductivity in La1.64Ca0.36Ga3O7.32 Melilite[END_REF] as a template. This had the double advantage of providing a good approximation of the underlying melilite framework and the oxide interstitial sites, without imposing symmetry constraints on the possible interstitial oxide ion orderings. We then performed a Rietveld refinement, allowing the atomic coordinates (52 independent atoms) to refine, along with the fractional occupancies of the interstitial oxide sites and three independent isotropic displacement parameters (one for each atom type La, Ga, O). Additionally, La4Ga2O9 was included as a secondary phase (note that this phase was only observed in the large-scale sample prepared for SPD-NPD combined refinement, at ~2wt%). This preliminary refinement produced a good fit to both data sets and revealed a distinct ordering of the interstitial oxides into chains (see Figure S4 for Rietveld fits and Table S2 for refined parameters). Having established a provisional model in this way, we sought to minimise the number of refined structural parameters in the model by performing a DFT optimisation of all atomic positions with lattice parameters fixed to experimental values, followed by a symmetry search. This transformed the structure into a new orthorhombic unit cell in space group Ima2, whose lattice parameters correspond to a √2a x √2a x 2c expansion of the tetragonal parent melilite unit cell. This model, which indexes both the SPD and NPD patterns correctly, contains 18 atoms in the asymmetric unit (vs. 50 in the provisional P1 model, Table S2) including 5 crystallographically inequivalent Ga atoms (vs. 12 in P1).

To confirm the space group assignment, we probed the number of crystallographically inequivalent Ga sites in the structure by using high-resolution solid-state NMR with 71 Ga (I = 3/2) as the active nucleus. [START_REF] King | A Space Group Assignment of ZrP2O7 Obtained by 31 P Solid State NMR[END_REF][START_REF] Martineau | Structure Resolution of Ba5Al3F19 and Investigation of Fluorine Ion Dynamics by Synchrotron Powder Diffraction, Variable-Temperature Solid-State NMR, and Quantum Computations[END_REF] As shown in Figure 2, the 71 Ga central-transition (CT) magic angle spinning (MAS) spectrum recorded at 20.0 T with a spinning frequency of 100 kHz exhibits a complex line shape revealing the presence of several overlapping 71 Ga resonances that are strongly affected by second-order quadrupolar effects. We therefore used a 2D double-quantum filtered (DQF) ST-MAS [START_REF] Kwak | Double-Quantum Filtered STMAS[END_REF] experiment to improve the spectral resolution by removing the second-order quadrupolar broadening in an indirect dimension. This allowed five distinct 71 Ga resonances to be resolved unambiguously, two of them (Ga1 and Ga2) with very close isotopic frequencies being discriminated by their different 71 Ga isotropic chemical shifts and second-order quadrupolar coupling parameters. The relative intensities of the resonances were obtained by fitting the quantitative CT MAS spectrum using the 71 Ga chemical shift and quadrupolar coupling parameters determined from the DQF-STMAS spectrum. This showed that the five resonances have relative intensities in the ratio 2 : 1.03 : 0.96 : 1.05 : 0.96, in full agreement with Ima2 symmetry for which the five crystallographically inequivalent Ga sites have a multiplicity ratio of 2:1:1:1:1. On this basis, we consider that symmetry lowering to P1 (leading to 12 inequivalent Ga sites of equal multiplicity) is not required. The final Rietveld refinement was carried out simultaneously against NPD and SPD data starting from the DFT-optimised Ima2 model. As shown in Figure 1(a,b), this produced good fits to both data sets (Rwp = 7.49 %, χ 2 = 3.06, which are nearly identical to the agreement factor χ 2 = 2.98 obtained in P1, as shown in Figure S4). The refined structure (final lattice parameters a = 11.4701(4) Å, b = 11.2674(4) Å, c = 10.4804(4) Å) is illustrated in Figure 3. To verify the model's reliability on the local scale, we conducted a simultaneous Rietveld-PDF fit against X-ray total scattering data (I15, Diamond Light Source, UK) using the Ima2 Rietveld model with fixed atomic coordinates. As shown in Figure 1(c,d), this provided a good fit to G(r) including over short distances of < 5Å (Rwp = 2.75 %). Finally, GIPAW [START_REF] Pickard | All-Electron Magnetic Response with Pseudopotentials: NMR Chemical Shifts[END_REF] DFT calculations of 71 Ga NMR parameters were performed for the refined Ima2 model. This produced calculated 71 Ga isotropic chemical shifts and quadrupolar coupling parameters that are consistent with the experimentally observed values (see Table 1).

Table 1:

Experimental and GIPAW-DFT calculated (italic) 71 Ga isotropic chemical shifts (CS ISO ), quadrupolar coupling constant (CQ), asymmetry parameter (ηQ) of the quadrupolar coupling tensor and relative intensities (I) for the Rietveld-refined Ima2 structural model. Calculated isotropic chemical shifts were deduced from calculated isotropic shieldings (σCS ISO ) using the relation CS ISO = -0.912 σCS ISO + 1519 (ppm). Experimental and GIPAW-calculated (italic) peak positions in the STMAS isotropic dimension (STMAS ISO ) for a spin-3/2 nucleus are given as STMAS ISO = CS ISO + 1 17

10 6 ν 0 2 ( 𝐶 𝑄 2 ) 2 (1 + η Q 2 3
), where 0 is the 71 Ga Larmor frequency (259.4 MHz). [START_REF] Gan | Isotropic NMR Spectra of Half-Integer Quadrupolar Nuclei Using Satellite Transitions and Magic-Angle Spinning[END_REF][START_REF] Amoureux | Unified Representation of MQMAS and STMAS NMR of Half-Integer Quadrupolar Nuclei[END_REF] The quadrupolar coupling parameters are defined as CQ =

eQ|V ZZ | ℎ and ηQ =| V YY _-V XX V ZZ
| where Q is the 71 Ga nuclear quadrupole moment, e is the electron charge, h is the Plank constant and Vii are the principal components of the EFG tensor at the Ga site (V.m 2 ) with |VXX| < |VYY| < |VZZ|. For each Ga site, the coordination polyhedron and its number of bridging oxygen atoms (BO) forming short Ga-O-Ga bonds are given. The refined structure of La2Ga3O7.5 retains the topology of the parent melilite LaAEGa3O7 (see Figure 3a,b). The interstitial oxide ions (Oint) are incorporated entirely into the [Ga3O7.5] layer, which is constructed from a fully-ordered corner-sharing network of GaO4 tetrahedra and GaO5 trigonal bipyramids in a 5:1 ratio. As in the parent structure, the [Ga3O7.5] framework forms pentagonal channels along the stacking axis. The La 3+ cations occupy the inter-layer space, where each is aligned to the centre of a pentagonal channel. This provides a suitable coordination environment for the interstitial oxide close to the centres of the channels (Figure 3c,d) where it is stabilised by significant bonding contributions from three neighbouring cations: a pair of short contacts to La 3+ ions from the adjacent layers (La2-Oint = 2.432( 5)Å and La3-Oint = 2.450(5)Å), and a short contact to the Ga 3+ at the apex of the pentagon to form a GaO5 unit (Ga5-Oint = 2.016(4)Å). These are supplemented by a pair of longer contacts to Ga 3+ (Ga1-Oint = 2.355(3)Å), and together these five contacts provide a calculated bond valence sum (BVS) of -1.81, with La-O and Ga-O contacts making a near-equal contribution. These bonding contributions between Oint and Ga 3+ cations are reflected in 71 Ga chemical shift trends: the Ga4 site (ISO~240 ppm) which is coordinated by 3 bridging (BO) and 1 non-bridging oxygen (NBO) atoms is shifted toward higher frequency with respect to the Ga2 (~198 ppm) and Ga3 sites (~214 ppm) which exhibit 4 Ga-O-Ga bridging bonds, as is expected since these 3 Ga sites are at remote distance from Oint. In contrast, the Ga5 site with 3 BO, 1NBO and an additional short contact to Oint exhibits an isotropic chemical shift of ~160 ppm typical of GaO5 units. The Ga1-Oint longer contact supplementing 3 bridging and 1 non-bridging bonds also leads to a significant shift toward lower frequency of the Ga1 site (δISO~187 ppm) relative to the Ga4 resonance. This mode of interstitial coordination is consistent with that observed in La1+xAE1-xGa3O7+x/2. [START_REF] Kuang | Interstitial Oxide Ion Conductivity in the Layered Tetrahedral Network Melilite Structure[END_REF][START_REF] Li | Interstitial Oxide Ion Order and Conductivity in La1.64Ca0.36Ga3O7.32 Melilite[END_REF][START_REF] Xu | La1+xBa1-xGa3O7+0.5x Oxide Ion Conductor: Cationic Size Effect on the Interstitial Oxide Ion Conductivity in Gallate Melilites[END_REF] All of the Oint sites are symmetrically equivalent. By calculating an O 2-BVS map from the refined structure, and searching for regions with BVS values between 1.8 -2.2 (shown as slices through the [Ga3O7.5] layers in Figure 3g,h) it can be seen that these regions coincide only with the fully-occupied Oint sites; the structure is distorted such that the unoccupied channels to not offer viable alternative coordination sites for Oint.

Site

In contrast to the x ≤ 0.65 melilites, La2Ga3O7.5 exhibits complete long-range ordering of Oint both within and between the layers. Within the layers, the occupied pentagonal rings form infinite vertex-linked chains which run parallel to the a-axis (Figure 3g,h). Each unit cell contains one chain per layer, so that the intra-layer distance between chains is equal to the b parameter (11.2674(4) Å). Between the layers, the chains are staggered to produce a 2-fold expansion of the c-axis relative to the tetragonal parent melilite. The staggering originates from interlayer ordering within the pentagonal channels: two distinct channel types can be identified and categorised as half-occupied (the La2/La3 channels, as shown in Figure 3e) or unoccupied (the La1 channels, Figure 3f) by Oint, with each unit cell containing four of each channel type. The occupancy of Oint in the La2/La3 channels alternates between layers such that each occupied Oint has a pair of empty neighbours. The alternation in Oint occupancy is coupled to an alternation in the inter-layer La-La distance: this provides the short La-La distance (4.882(1) Å) required for the local coordination of Oint, and an alternating long La-La distance (5.600(1) Å) which traverses the empty adjacent sites. In contrast, the unoccupied channels (La1, Figure 3f), which have no such local coordination requirements, exhibit a regular La-La distance of 5.2404(3) Å (= c/2). The chain-like Oint ordering in La2Ga3O7.5 is consistent with the local occupancy rule proposed by Li et al. for La1.64Ca0.36Ga3O7.32, i.e. that the centres of the five pentagonal rings adjacent to Oint cannot themselves be occupied, due to the local strain produced by framework relaxation. [START_REF] Li | Interstitial Oxide Ion Order and Conductivity in La1.64Ca0.36Ga3O7.32 Melilite[END_REF] However, chain-like ordering is not the only way to satisfy this rule: two-dimensional ordering (the square arrangement of Oint proposed by Li et al.) [START_REF] Li | Interstitial Oxide Ion Order and Conductivity in La1.64Ca0.36Ga3O7.32 Melilite[END_REF] provides an alternative model which is intuitive because it maximises the Oint-Oint separation. The derivation of the chain-and square-orderings from this local occupancy rule is illustrated in Figure S6. To determine the relative stabilities of these two orderings in La2Ga3O7.5, we generated a 1x2x1 supercell for each structure and relaxed them in P1 symmetry by DFT (see Figure 4b,c). This confirmed that chain-like ordering produces the most energetically favourable structure, with a difference of 0.44eV between the two models. The difference between the stabilities of the two ordering models can be rationalized by representing the two dimensional [Ga3O7] layer as a pair of interpenetrated hexagonal nets ("Cairo pentagonal tiling" or "MacMahon's Net") [START_REF] O'keeffe | Plane Nets in Crystal Chemistry[END_REF] as shown in Figure 4a. This means that the [Ga3O7] layer can be considered as a hexagonal tiling, with a single tile whose vertices coincide with the Ga1 atom positions. The tile contains four edgesharing pentagons as a motif (thus generating the orthogonal hexagonal net when tiled, as shown in Figure 4d). The inclusion of one interstitial oxide ion in the motif, at the center of one of the four pentagons generates the correct stoichiometry for the layer, [Ga3O7.5]. Finally, two edge-lengths d and d' are defined, which correspond respectively to the Oint-occupied and Oint-unoccupied halves of the hexagonal tile (Figure 4d). Using this tile, it is possible to generate the different orderings by antiparallel tiling (chain ordering) or parallel tiling (square ordering), as shown in Figure 4d. Inspection of these tilings reveals that the two different ordered structures do not have the same freedom to distort locally. Specifically, parallel tiling (square ordering) is subject to a geometrical constraint that requires d = d'. This places a constraint on the local geometry of the [Ga3O7.5] framework, as confirmed by inspection of the d and d' distances (Ga1 -Ga1) in the DFTrelaxed P1 supercell, which are near-equal (6.337 and 6.338Å, Figure 4c). In contrast, antiparallel tiling (chain ordering) is subject only to the general condition that the hexagonal tile have two parallel edges of equal length. Consequently, a distortion with d ≠ d' is permitted, and inspection of these distances relaxed by DFT in P1 confirms the occurrence of a local framework relaxation of this type (6.174 and 6.481Å, Figure 4b). In summary, the chain-like ordering of Oint within the [Ga3O7.5] layers of La2Ga3O7.5 results from the underlying hexagonal topology of the melilite framework, and allows the framework to have the maximum flexibility in response to the local strain produced by incorporation of these interstitial ions. In contrast, the two dimensional square ordering proposed on a local scale for La1.64Ca0.36Ga3O7.32 [START_REF] Li | Interstitial Oxide Ion Order and Conductivity in La1.64Ca0.36Ga3O7.32 Melilite[END_REF] would impose additional geometrical requirements which prevents the framework from relaxing fully when interstitial defects are introduced. These observations also cast new light on the complex partially-ordered structure of La1.64Ca0.36Ga3O7.32, as shown in Figure S7 which suggests that chain-like ordering is incipient in this compound.

Stability and structure at high temperature. We collected PXRD data at elevated temperatures to evaluate the thermal stability of La2Ga3O7.5. The compound was found to be stable below 830°C: above this temperature it begins to decompose into a mixture of LaGaO3 and Ga2O3 (Figure 5a). La2Ga3O7.5 remains metrically orthorhombic throughout its stability range, in contrast to the orthorhombic melilites La1.64Ca0.36Ga3O7.32 and La1.64Sr0.36Ga3O7.32, which undergo a reversible phase transition to a disordered tetragonal polymorph on heating. [START_REF] Li | Interstitial Oxide Ion Order and Conductivity in La1.64Ca0.36Ga3O7.32 Melilite[END_REF][START_REF] Thomas | Phase Stability Control of Interstitial Oxide Ion Conductivity in the La1+xSr1-xGa3O7+x/2 Melilite Family[END_REF] An incipient phase transition of this type is evident in La2Ga3O7.5 as indicated by the trend in its basal lattice parameters a and b: extrapolation of the trend in a/b ratio by fitting to a power law of form 𝑓(𝑇) = 𝐴 -𝐵. 𝑇 -(𝑇 𝐶 -𝑇 ) -𝐶 indicates a hypothetical phase transition temperature of ~890°C, which is higher than the decomposition temperature. However, the onset of this transition occurs at much lower temperatures, as suggested by the deviation from linearity in the range 610 -710°C (Figure 5c). The thermal expansion coefficients calculated from the linear region (T ≤ 610°C), shown in Table S4, are consistent with those reported for La1+xSr1-xGa3O7+x/2. 12 Conductivity. Figure 6a shows the complex impedance plot of a La2Ga3O7.5 disk (obtained by polishing a bead) at 400 °C in an oxygen atmosphere, which comprises one large semicircular arc and an inclined line. The large semicircular component of the high frequency arc displays capacitance ~ 4-6 pF cm -1 , and is ascribed to the bulk response, which can be simply modelled with a parallel resistor (Rb) and capacitor (Cb). The intercept of this arc at low frequency was extracted as bulk resistivity Rb and the calculated capacitance Cb is ~ 4 pF cm -1 according to the equation 2fmax RC = 1 (fmax is the frequency at maximum imaginary impedance Z''max). This value is close to that calculated using the equation M''max = 0/2C (M''max is the maximum imaginary modulus), confirming that the large semicircular arc is from the single bulk response. The inclined line in the low frequency region has capacitances ranging from 10 -11 F cm -1 to 10 - 7 F cm -1 (inset of Figure 6a), suggesting heavily overlapping grain boundary and electrode response arcs. At high temperature, the bulk response arc gradually disappeared and the grain boundary and electrode response arcs became well resolved (see 700°C data, Figure 6b). The large capacitances 10 -7 ~ 10 -5 F cm -1 for the electrode response arc suggest the occurrence of ionic conduction in La2Ga3O7.5.

Figure 6c shows the Arrhenius plot of bulk conductivities of La2Ga3O7.5. In the range 300-800 °C, the bulk conductivity for La2Ga3O7.5 increases from 2.9×10 -8 -0.01 S cm -1 , which is three orders of magnitude higher than the parent compound LaSrGa3O7. However, La2Ga3O7.5 is less conductive than La1.54Sr0.46Ga3O7.27. Below 600°C, the conductivity of La2Ga3O7.5 is three orders of magnitude lower than that for La1.54Sr0.46Ga3O7.27, and has an activation energy of 1.21(3) eV. At high temperatures above 600 °C, the difference between conductivities of La2Ga3O7.5 and La1.54Sr0.46Ga3O7.27 reduces to approximately one order of magnitude, but the apparent activation energy of La2Ga3O7.5 increases to 1.62(2) eV. Although the VT-PXRD data shows the incipient phase transition from orthorhombic to tetragonal starts at temperatures above 610 °C, no typical abrupt jump in the conductivity was observed over the measured temperature region before the sample starts to decompose. We ascribe this behavior to a gradual thermally-induced disordering of oxide interstitials -a process that does not reach completion -which enhances the oxide ion conductivity and increases superficially the measured activation energy in the high temperature region. A lower activation energy region would be expected to appear in a higher temperature region for the hypothetical fully-disordered tetragonal phase if the material survived at such high temperatures: such behavior is exhibited in the well-defined phase transition region of La1.64Ca0.36Ga3O7.32 where a steeper gradient in the Arrhenius plot at ~600°C links the low-and high-temperature plateaus associated with the orthorhombic and tetragonal phases respectively. 18

CONCLUSIONS

A new ternary melilite compound La2Ga3O7.5 can be synthesized by direct crystallization of an under-cooled melt, showing that it is possible to exceed the long-established limit of interstitial oxide concentration in the melilite structure by selecting a suitable non-equilibrium synthesis method. La2Ga3O7.5 represents the end-member (x = 1) of the La1+xAE1-xGa3O7+x/2 interstitial oxide ion conductor family, and is the first example of a ternary melilite compound. It adopts a unique anion-ordered superstructure with chain-like ordering of interstitial oxide ions, which is driven by the interplay between the local coordination requirements of the interstitial oxides and the underlying hexagonal topology of the [Ga3O7] host framework. La2Ga3O7.5 has a conductivity of 0.01 Scm -1 at 800°C and is stable on heating up to 830°C, putting it in the range of interest for solid oxide fuel cell devices. Despite the retention of long-range structural ordering up to its decomposition temperature, this is a reasonably high conductivity, three orders of magnitude above than that of the corresponding x = 0 end member LaSrGa3O7 and only one order of magnitude below that of La1.54Sr0.46Ga3O7.27. These observations will open the door to extended ranges of La1+xAE1-xGa3O7+x/2 melilite conductors with the potential for increasing conductivities in the Ca, Sr, Ba series, and may have a broader impact on related interstitial oxide conductor families such as the langasites and apatites, where the exploration of extended compositional ranges could lead to the discovery of new high-performance materials.

Figure 1 .

 1 Figure 1. Structural refinements of La2Ga3O7.5 in Ima2. (a,b) Combined Rietveld refinement against SPD data from 11BM, and NPD data from D2B, with the fits to the high-Q regions inset. (c,d) Simultaneous Rietveld and PDF fit against SPD data from I15-1. Magenta ticks = La2Ga3O7.5, cyan ticks = La4Ga2O9 (~2wt%). Black points = yobs, red line = ycalc, blue line = yobs -ycalc.

Figure 2 :

 2 Figure 2: (a) Experimental (black) and reconstructed (red) 71 Ga quantitative CT MAS spectra of La2Ga3O7.5 recorded at 20.0 T with a spinning frequency of 100 kHz. The asterisks indicate spinning sidebands. The five individual 71 Ga resonances are shown below. (b) Experimental (colored contours) and reconstructed (gray contours) DQF STMAS spectra obtained at 20.0 T with a spinning frequency of 100 kHz evidencing the five distinct 71 Ga resonances.

Figure 3 .

 3 Figure 3. Refined crystal structure of La2Ga3O7.5 in Ima2. (a,b) Unit cell projections along the aand c-axes (blue = La 3+ , green = Ga 3+ , red = O 2-). (c,d) Local coordination environment of the interstitial oxide at the centres of the pentagonal rings, labelled with distances (in Å) for the nearest O-Ga and O-La contacts. (e,f) The two distinct channel types labelled with their inter-layer La-La distances (in Å), showing a periodic short/long alternation in the interstitial-hosting channels (the La2/La3 channels) which is commensurate with the interstitial oxide ordering, in contrast to the regular La-La spacing in the "empty" channels (the La1 channels). (g,h) calculated BVS maps for O 2-, sliced

Figure 4 .

 4 Figure 4. DFT comparison of experimentally-observed chain-like oxide ion ordering and hypothetical square ordering. (a) The melilite [Ga3O7] framework projected along 00l, viewed as two interpenetrating hexagonal nets. Green = Ga, red = O. (b, c) DFT-relaxed 1x2x1 P1 supercells of La2Ga3O7.5 with interstitial oxide arranged into chains or squares, overlaid with a fragment of the underlying hexagonal net. Green = GaO4, magenta = GaO5, blue = La. (d) Derivation of the chain-and square orderings from antiparallel and parallel tilings of a single hexagonal tile. Note that parallel tiling (square ordering) imposes a special geometric constraint on the tile (d = d'), corresponding to a loss of [Ga3O7.5] framework flexibility when compared to the antiparallel chain-ordered case.

Figure 5 :

 5 Figure 5: Structural trends in La2Ga3O7.5 on heating. (a) PXRD patterns collected on heating from 30 -850°C showing the onset of decomposition to a perovskite phase (peaks marked with *) above 830°C. (b) Relative trends in the lattice parameters, plotted as y/y0-1 (where y = a, b, c or V and y0 are their corresponding room temperature values). (c) Relative trend in the orthorhombicity of the unit cell, plotted as a/b-1, and fitted to a power law of form 𝑓(𝑇) = 𝐴 -𝐵. 𝑇 -(𝑇 𝐶 -𝑇 ) -𝐶 (red line; see TableS3for fitted values of A, B, C and the critical temperature TC). The black dashed line is a guide to the eye that follows the initial linear trend from 30°C. The gray shaded areas in (b) and (c) indicate the region above the melilite decomposition temperature.

Figure 6 .

 6 Figure 6. Complex impedance plots at (a) 400 °C and 700 °C of the La2Ga3O7.5 ceramic. Rb and Rb+Rgb denote the bulk and total resistivities, respectively. The inset in each part shows capacitance versus frequency plot. (c) Arrhenius plot of the bulk conductivities for La2Ga3O7.5 (red points) in comparison with those of LaSrGa3O7 (black points) and La1.54Sr0.46Ga3O7.25 (black line, data from reference 14).

ACKNOWLEDGMENT

Financial support from the CNRS (IR-RMN-THC Fr3050), the State Scholarship Fund of China Scholarship Council (201704910448), the Natural Science Foundation of China (No. 21622101), and Guangxi Natural Science Foundation (No. 2019GXNSFGA245006) are gratefully acknowledged. We thank the beamline staff at 11BM (APS, Argonne National Laboratory, USA) for collecting the SPD data via remote access, Dr. E. Suard (ILL, Grenoble, France) for performing the neutron diffraction experiment via EASY access, and Dr. E. Veron (CNRS, CEMHTI UPR3079, Orléans) for performing the in-house SEM-EDX measurements.

ASSOCIATED CONTENT

Supporting Information. Details of attempts to synthesize La2Ga3O7.5 by low temperature glass-crystallization; Figures S1 -S7, Tables S1 -S4; crystallographic information file. This material is available free of charge via the Internet at http://pubs.acs.org.

AUTHOR INFORMATION

Corresponding Author * michael.pitcher@cnrs-orleans.fr * mathieu.allix@cnrs-orleans.fr


Yurena Arjona
email: yurenaaf@gmail.com

Javier Fernández-López

| Miguel Navascués

| Nadir Alvarez

Manuel Nogales

Pablo Vargas

Ceridwen Fraser

Manuel Noga

Linking seascape with landscape genetics: Oceanic currents favour colonization across the Galápagos Islands by a coastal plant

Aim: Coastal plants are terrestrial organisms for which ocean surface currents often act as long-distance dispersal vectors (thalassochorous species) favouring broad distributions and connecting distant populations. However, few studies have statistically assessed the role of currents in modulating gene flow and species distributions of terrestrial organisms. Here we evaluate the hypothesis that some thalassochorous plants exhibit population connectivity, presumably due to effective seed dispersal driven by sea currents.

Location: Galápagos Islands (Ecuador).

Taxon: Salt bush (Cryptocarpus pyriformis Kunth), a Galápagos native and locally widespread coastal angiosperm.

Methods: Using 1806 SNPs obtained by ddRADseq, we evaluated the genetic structure and differentiation of the Galápagos salt bush. To assess the role of sea currents in modulating inter-population gene flow, four explicit hypotheses were tested using reciprocal causal modelling and spatial eigenvector analysis: (a) isolation by sea resistance, considering that only sea dispersal is possible; (b) isolation by sea and inland resistance, considering that inland dispersal is also possible; (c) isolation by barrier, considering the sea as an obstacle to seed dispersal; and (d) isolation by geographical distance.

Results: Low differentiation and little genetic structure were detected among populations of C. pyriformis. Pairwise genetic distances between populations from different islands were significantly correlated with cost distances calculated from sea-current direction and speed. Nonetheless, inland dispersal also accounted for some gene flow within each island.

Main conclusion:

Extensive and frequent seed dispersal by sea has apparently favoured strong inter-island genetic connectivity within Galápagos. A combination of methods developed for terrestrial and marine domains (landscape and seascape genetics) aids in understanding how landscape features modulate gene flow of coastal plant species, as these terrestrial organisms are highly dependent on the sea for seed dispersal.

| INTRODUC TI ON

Oceans cover around 71% of the Earth surface, isolating landmasses and hence the terrestrial organisms that occupy them. At the same time, sea currents act as an essential vector for propagules of land organisms, transporting them across long distances, connecting inhabited landmasses and allowing colonization of oceanic islands that emerge lifeless from the seafloor [START_REF] Abe | Colonization of Nishino-Shima island by plants and arthropods 31 years after eruption[END_REF][START_REF] Magnússon | Plant colonization, succession and ecosystem development on Surtsey with reference to neighbouring islands[END_REF][START_REF] Ridley | The dispersal of plants throughout the world[END_REF]. Plant species specialized for sea dispersal are referred to as thalassochorous.

They are able to develop floating diaspores that survive a long time in seawater, which favours their arrival and germination in new territories [START_REF] Van Der Pijl | Principles of dispersal in higher plants[END_REF][START_REF] Ridley | The dispersal of plants throughout the world[END_REF].

Thalassochorous plants usually have wide distributions [START_REF] Miryeganeh | Long-distance dispersal by sea-drifted seeds has maintained the global distribution of Ipomoea pes-caprae subsp. brasiliensis (Convolvulaceae)[END_REF], which points to the efficiency of the thalassochorous dispersal syndrome in long-distance dispersal (hereafter LDD). Coastal plants, due to their proximity to the sea, are salt tolerant and many of them develop diaspores with corky or fibrous tissues, air chambers and a diverse array of structures that help them float for long periods over the sea. This set of traits defines the thalassochorous dispersal syndrome (Heleno & Vargas, 2015;[START_REF] Van Der Pijl | Principles of dispersal in higher plants[END_REF][START_REF] Ridley | The dispersal of plants throughout the world[END_REF]. [START_REF] Darwin | On the origin of species by means of natural selection, or the preservation of favoured races in the struggle for life[END_REF] experimentally demonstrated the ability of certain plant species to float for a long time and germinate after immersion in salt water, which has been confirmed by subsequent studies reporting a wide variety of coastal species that are able to float for months in salt water and germinate afterwards [START_REF] Esteves | On the limited potential of Azorean fleshy fruits for oceanic dispersal[END_REF][START_REF] Guja | Buoyancy, salt tolerance and germination of coastal seeds: Implications for oceanic hydrochorous dispersal[END_REF][START_REF] Nakanishi | Dispersal ecology of the maritime plants in the Ryukyu Islands[END_REF]. Such high dispersal ability and salt tolerance of thalassochorous plants have favoured the early colonization of oceanic islands [START_REF] Abe | Colonization of Nishino-Shima island by plants and arthropods 31 years after eruption[END_REF][START_REF] Magnússon | Plant colonization, succession and ecosystem development on Surtsey with reference to neighbouring islands[END_REF], the recolonization of disturbed areas and genetic connectivity between distant populations by the exchange of diaspores [START_REF] Cerón-Souza | Contrasting demographic history and gene flow patterns of two mangrove species on either side of the Central American Isthmus[END_REF][START_REF] Takayama | Gene flow and population subdivision in a pantropical plant with sea-drifted seeds Hibiscus tiliaceus and its allied species: Evidence from microsatellite analyses[END_REF]. In principle, gene flow between sea-dispersed plant populations should be highly influenced by sea current properties, such as directionality and speed, which should be major drivers in shaping their population genetic structure [START_REF] Gallaher | Navigating the 'broad freeway': Ocean currents and inland isolation drive diversification in the Pandanus tectorius complex (Pandanaceae)[END_REF].

Understanding the effect of landscape features on gene flow between populations of terrestrial organisms is one of the objectives of landscape genetics [START_REF] Manel | Ten years of landscape genetics[END_REF]. In contrast to terrestrial organisms, marine organism gene flow is highly influenced by seascape features. As a consequence, seascape genetics was developed as a discipline closely related to landscape genetics to cope with the particular biology of many marine organisms-that is, wide distributions, high dispersibility and little genetic structureand physical properties of the ocean [START_REF] Riginos | Seascape genetics: Populations, individuals, and genes marooned and adrift[END_REF][START_REF] Selkoe | A decade of seascape genetics: Contributions to basic and applied marine connectivity[END_REF]. Coastal plant species are in an intermediate situation, being terrestrial organisms that are mainly sea dispersed but capable of dispersal over land. Although many phylogeographical studies suggest that ocean currents are responsible for the genetic structuring patterns found (e.g. [START_REF] Gallaher | Navigating the 'broad freeway': Ocean currents and inland isolation drive diversification in the Pandanus tectorius complex (Pandanaceae)[END_REF][START_REF] Miryeganeh | Long-distance dispersal by sea-drifted seeds has maintained the global distribution of Ipomoea pes-caprae subsp. brasiliensis (Convolvulaceae)[END_REF][START_REF] Takayama | Gene flow and population subdivision in a pantropical plant with sea-drifted seeds Hibiscus tiliaceus and its allied species: Evidence from microsatellite analyses[END_REF], few studies have specifically addressed the extent to which direction and speed of sea currents modulate gene flow among coastal plant populations (e.g. see [START_REF] Hodel | Terrestrial species adapted to sea dispersal: Differences in propagule dispersal of two Caribbean mangroves[END_REF]. A combination of methods employed in landscape and seascape genetics is the best option to address this issue.

Cryptocarpus pyriformis Kunth (Nyctaginaceae) -commonly known as salt bush or monte salado -is distributed along the coasts of Galápagos, mainland Ecuador and Peru [START_REF] Wiggins | Flora of the Galápagos Islands[END_REF]). The Galápagos archipelago, located in the Pacific Ocean about 1,000 km west of mainland Ecuador, hosts many populations of salt bush, occurring on all 13 largest islands [START_REF] Jaramillo Díaz | CDF checklist of galapagos flowering plants[END_REF]. The small flowers of C. pyriformis are mainly pollinated by insects, with ants being the most frequent visitor of the plants in the Galápagos [START_REF] Traveset | Invaders of pollination networks in the Galápagos Islands : Emergence of novel communities[END_REF]. Its fruits are small, pear-shaped and have a corky consistency that helps them to float on salt water (Y. Arjona and P. Vargas, pers. obs.). Galápagos finches feed on salt bush fruits and intact seeds have been found in lava lizard droppings, although no conclusive results have been reported regarding effective dispersal by vertebrates [START_REF] Heleno | Frugivory and seed dispersal in the Galápagos: What is the state of the art?[END_REF][START_REF] Heleno | Seed dispersal networks in the Galápagos and the consequences of alien plant invasions[END_REF]. There is thus evidence of a dispersal potential over land as well as dispersal by sea in this species.

In this study, we evaluate the hypothesis that coastal thalassochorous species readily disperse by sea, which should be reflected by evidence of gene flow corresponding to the direction and speed of sea currents. With this aim we chose the widespread salt bush (Cryptocarpus pyriformis) as a case study and evaluated its dispersal across the Galápagos archipelago. We formulated four explicit hypotheses that consider different scenarios that vary in importance for inland and sea dispersal: (a) dispersal is largely mediated by sea and modulated by direction and speed of sea currents, whereas inland dispersal does not contribute to inter-population connectivity (isolation by sea resistance); (b) sea dispersal is the main driver of inter-island connectivity, and inland dispersal also contributes to genetic differentiation (isolation by inland and sea resistance); (c) the sea acts as a barrier to dispersal, isolating populations from different islands (isolation by barrier); and (d) isolation by geographical distance, with no effect of the context (island or sea; null hypothesis).

Publicly available data on ocean circulation and genetic data-from hundreds of single nucleotide polymorphisms (SNPs) distributed across the whole genome-were analysed under a landscape genetic framework, to evaluate the role of sea current speed and direction in inter-island genetic differentiation of the salt bush across the Galápagos.

K E Y W O R D S

directional dispersal, isolation by resistance, oceanic islands, population genetics, salt bush, thalassochory

| MATERIAL S AND ME THODS

| Sampling, DNA extraction, ddRAD library preparation and SNP calling

A total of 231 samples were collected from 12 of the 13 largest islands of the Galápagos (Baltra was excluded because of proximity to Santa Cruz and its recent land bridge), where we selected from one to five populations per island based on island size. Additionally, nine individuals were sampled from four populations along the coast of mainland Ecuador, to assess how divergent the Galápagos populations are from distant continental populations (Figure 1 and Table S1.1). DNA was extracted from silica-gel-dried leaf samples using DNeasy Plant Mini Kit (Qiagen, Inc., Valencia, CA, USA) and cleaned following the phenol-chloroform protocol to eliminate salts and other compounds that may interfere in further steps of the process. Double-digest restriction site associated DNA (ddRAD) libraries were prepared following the protocol described by [START_REF] Mastretta-Yanes | Restriction site-associated DNA sequencing, genotyping error estimation and de novo assembly optimization for population genetic inference[END_REF], including 24 replicates. The ddRAD libraries were single-end sequenced on two lanes of Illumina HiSeq 2,500 at the Lausanne Genomic Technologies Facility.

Reads were demultiplexed and de novo assembled using 'ipyrad' v. 0.7.22 [START_REF] Eaton | PyRAD: Assembly of de novo RADseq loci for phylogenetic analyses[END_REF]. 'Trimmomatic' v. 0.36 (Bolger, Lohse, & Usadel, 2014) was used to remove the adapters and quality-trim the reads before the assembly. The 'ipyrad' clustering threshold and minimum depth parameters were set to 0.93 and 10, respectively, as they rendered a good balance between the error rate and the number of SNPs recovered (Figure S2.1). The resultant SNP matrix was filtered using 'vcftools' v. 0.1.15 [START_REF] Danecek | The variant call format and VCFtools[END_REF]; all sites with more than 70% and all individuals with more than 85% F I G U R E 1 Galápagos archipelago and mean values of sea current speed and direction for the whole period analysed . Arrows indicate direction, their length and width are proportional to the speed of the sea currents. Location of the Galápagos archipelago is shown in the map (top right, mainland Ecuador in dark grey). Red dots indicate the position of Cryptocarpus pyriformis sample populations. CO: continent; ES: Española; SC: San Cristóbal; FL: Floreana; SF: Santa Fe; SX: Santa Cruz; PZ: Pinzón; SA: Santiago; IS: Isabela; FE: Fernandina; GE: Genovesa; MA: Marchena; PT: Pinta of missing data were removed, and only biallelic sites were kept.

Additionally, custom scripts were employed to further modify 'ipyrad' output files to fit the input requirements of the software employed when necessary. See Appendix S2 for further details about ddRAD library preparation and bioinformatic processing of the reads.

| Genetic diversity and population genetic structure

Observed and expected heterozygosities (H O , H E ) and the inbreeding coefficient (F IS ) were calculated across loci per population using 'GenAlEx' [START_REF] Peakall | GenAlEx 6: Genetic analysis in Excel. Population genetic software for teaching and research[END_REF], 2012). A Pearson's correlation was performed between heterozygosity and island size to evaluate a possible variation in genetic diversity based on different carrying capacities due to island size. Jost's D differentiation index [START_REF] Jost | GST and its relatives do not measure differentiation[END_REF], multilocus fixation index G ST [START_REF] Nei | Estimation of fixation indices and gene diversities[END_REF], Weir and Cockerham's F ST estimator-θ [START_REF] Weir | Estimating F-Statistics for the analysis of population structure[END_REF])-across loci and population, and pairwise θ between populations, were estimated with the package 'diveRsity' [START_REF] Keenan | DiveRsity: An R package for the estimation and exploration of population genetics parameters and their associated errors[END_REF] for the r software (R Development Core Team, 2017). Analyses of molecular variance (AMOVA; [START_REF] Excoffier | Analysis of molecular variance inferred from metric distances among DNA haplotypes: Application to human mitochondrial DNA restriction data[END_REF] were performed to assess the proportion of variation between islands and between populations within islands using the r package 'poppr' [START_REF] Kamvar | Poppr : An R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction[END_REF].

To evaluate the degree to which genetic differentiation and genetic diversity increase when including very distant populations, we repeated the analysis including continental samples.

A principal component analysis (PCA) and a discriminant analysis of principal components (DAPC) were employed to explore the genetic variation among sampled individuals and the existence of genetically distinct groups of individuals. The groups assessed by DAPC can either be given by the k-means algorithm or defined by the user. We ran DAPC using both options, considering the 12 islands and the continent as groups defined by the user (K = 13). The number of principal components (PCs) retained in the analyses was selected based on the optimization of the α-score to avoid over-fitting the results. To perform PCA and DAPC we used the r package 'adegenet' [START_REF] Jombart | Discriminant analysis of principal components: A new method for the analysis of genetically structured populations[END_REF]. In addition, we used ' STRUCTURE' v. 2.3.4 (Falush, Stephens, & Pritchard, 2003;[START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF] to assess population genetic structure. We ran the admixture model assuming correlated allele frequencies, with a lambda previously estimated as recommended by the authors, from K = 1 to K = 13, 10 times with 2,000,000 MCMC iterations and 200,000 iterations burn-in. To assess which value of K was statistically preferred we used two methods: (a) ΔK [START_REF] Evanno | Detecting the number of clusters of individuals using the software STRUCTURE: A simulation study[END_REF]; and (b) MedMeaK, MaxMeaK, MedMedK and MaxMedK [START_REF] Puechmaille | The program STRUCTURE does not reliably recover the correct population structure when sampling is uneven: Subsampling and new estimators alleviate the problem[END_REF]. These latter estimators were calculated using different thresholds (0.5, 0.6, 0.7 and 0.8), considering all the 23 Galápagos populations and the continent. The graphical output from the 'STRUCTURE' results was created using 'DISTRUCT' v.1.1 [START_REF] Rosenberg | DISTRUCT: A program for the graphical display of population structure[END_REF].

| Resistance due to ocean currents within Galápagos

Resistance to seed dispersal due to ocean currents was estimated using publicly available data from the U.S. Naval Research Laboratory, specifically the HYCOM + NCODA Global 1/12º Analysis (GLBu0.08) dataset (https://www.hycom.org/data/glbu0 pt08). Data for the working area (-92.5º and -88.25º east; -2º and1.5º north) were obtained via functions from the r package 'HMMoce' [START_REF] Braun | HMMoce: An R package for improved geolocation of archival-tagged fishes using a hidden Markov method[END_REF]; https://github.com/cran/HMMoc e/ blob/maste r/R/get.hycom.R). These data consisted of the eastward and northward components of seawater velocity (m/s) at the ocean surface for 4,771 days (measures taken every two days during the period from 2 October 1992 to 20 November 2018, see Figure 1 for a summary of these data). This extensive period of time was selected to include annual and seasonal variations in climatic conditions and ocean circulation, and because Cryptocarpus fruits are produced throughout the year.

Anisotropic cost-distance matrices between sampling localities were estimated based on the maximum connectivity within the period covered by the data. Firstly, for each day, raw data were converted into a raster layer from which the speed and direction of sea currents were calculated. Then, the cost to passively move between adjacent cells of the raster was calculated using the function flow.dispersion from the 'rWind' package [START_REF] Fernández-López | rWind: Download, edit and include wind data in ecological and evolutionary analysis[END_REF], modified to cope with missing data values corresponding to landmasses. This function implements the algorithm proposed by [START_REF] Muñoz | Wind as a long-distance dispersal vehicle in the southern hemisphere[END_REF], in which the cost is directly proportional to the difference between the direction of the movement trajectory towards the target cell and the direction of the sea current, and inversely proportional to sea current speed. For passive dispersal, movements against sea current flow are not permitted-that is, those for which the difference between their bearing and sea current direction is higher than 90º (see [START_REF] Fernández-López | rWind: Download, edit and include wind data in ecological and evolutionary analysis[END_REF]. After their computation, transition costs among sea surface cells ranged between 0.1 and 32,004.9, with a mean of 310.7. Based on these costs, two resistance layers were calculated:

(a) considering only sea dispersal is possible (infinite cost assigned to inland dispersal); and (b) that both inland and sea dispersal occur (a resistance of 0.1 assigned to all land cells-this is the minimum cost to move between adjacent cells assigned by the function, without considering sea current speed). A third resistance layer was built to account for the hypothesis of isolation due to barriers, where infinite cost was assigned to sea dispersal and a cost of 0.1 to the transition between adjacent land cells. With the transition costs between all adjacent cells, least-cost paths between populations were calculated using the 'gdistance' R package [START_REF] Van Etten | R package gdistance: Distances and routes on geographical grids[END_REF]. From the inverse of these costs we obtained the connectivity between pairs of populations. This procedure yielded 4,771 connectivity matrices, one for each day in the dataset, summarized into a single matrix per hypothesis by taking the maximum value per population pair and direction and converted them again into cost-distance values. To make this asymmetric matrix symmetric (i.e. to have costs from A to B and from B to A summarized into a single distance A-B), we calculated the mean value between the two connectivity values per population pair. Independently, additional cost distance matrices were calculated under each of the two resistance hypotheses for the rainy season (from December to May) and the dry season (from April to November) separately. Due to the proximity of Española's populations ES1 and ES3, the estimated connectivity between them was infinite. To avoid this extreme value affecting the results, we removed population ES3 for further analyses.

| Isolation by resistance, barrier or geographical distance

We proposed four alternative models that modulate gene flow be-

| Reciprocal causal modelling

Pairwise genetic distances were calculated as F ST /(1 -F ST ) using 'Genepop' 4.0 [START_REF] Rousset | GENEPOP'007: A complete re-implementation of the GENEPOP software for Windows and Linux[END_REF] as recommended by [START_REF] Rousset | Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance[END_REF].

Geographical distances were calculated using the r package 'raster' [START_REF] Hijmans | raster: Geographic data analysis and modeling[END_REF] as the natural logarithm of the distance between two points in the WGS ellipsoid to linearize the relationship with genetic distances. Cost distances mediated by sea currents for the two IBR and the IBB models were obtained as explained in the previous section, and were also log-transformed. The correlations between distance matrices were calculated using Mantel test [START_REF] Mantel | The detection of disease clustering and a generalized regression approach[END_REF] and partial Mantel test [START_REF] Smouse | Multiple regression and correlation extensions of the Mantel test of matrix correspondence[END_REF]) using the 'vegan' R package [START_REF] Oksanen | vegan: Community Ecology package[END_REF]. Partial Mantel statistics (r M ) were then compared under reciprocal causal modelling, which reduces the high risk (reported for simple Mantel tests) of finding spurious correlations when comparing highly correlated alternative hypotheses [START_REF] Cushman | Re-evaluating causal modeling with Mantel tests in landscape genetics[END_REF]. For a pair of hypotheses (H1 and H2), the relative support of H1 is given by the difference between r M obtained from two partial Mantel tests:

(a) between genetic distances and H1, controlling for the effects of H2; and (b) between genetic distances and H2, controlling for the effects of H1 (i.e. H1|H2-H2|H1). If H1 is correct this difference is positive, otherwise H2 is supported. After repeating this procedure for all hypothesis pairs, a full matrix is obtained in which a hypothesis is fully supported when all values in the column are positive and the values in the associated row are negative (see Appendix S3 for a more detailed description of the method). This analysis was repeated considering the rainy and dry season separately, to evaluate the possibility of seasonal differences in connectivity that may be masked when analysing the whole dataset. In addition, to evaluate the effect of inter-annual variation due to the El Niño phenomenon on population connectivity, we estimated the correlation between El Niño Southern Oscillation Index (SOI) and the Mantel statistic for the correlation between genetic and cost distances calculated for each month of the data period (Appendix S4).

| Spatial eigenfunction analysis

Spatial eigenfunction analysis was performed as an alternative method to assess the relative contribution of geography and sea currents to the spatial genetic variation. The main difference with reciprocal causal modelling is the use of raw genetic data and redundancy analysis, instead of the pairwise genetic differences and Mantel tests that have lower statistical power [START_REF] Legendre | Should the Mantel test be used in spatial analysis?[END_REF]. The geographical distance matrix and the three asymmetric cost distance matrices-from IBR1, IBR2 and IBB hypotheses-were decomposed into eigenvector maps: (a) distance-based Moran's eigenvector maps (dbMEM) were obtained by decomposing pairwise Euclidean distances between populations (IBD hypothesis);

(b) Moran's eigenvector maps (MEM) were obtained from a spatial weighting matrix built from a neighbourhood matrix where only populations from the same island were considered neighbours, weighted by the pairwise cost distances between neighbours (IBB hypothesis) [START_REF] Dray | Spatial modelling: A comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM)[END_REF]; and (c) asymmetric eigenvector maps (AEM) as the result of decomposing the directional spatial process of sea flow (IBR1 and IBR2 hypotheses). Directional connection diagrams were obtained from the sea-cost matrices built under each resistance hypothesis. For each population pair, only the direction that resulted in the highest connectivity between the populations was selected, setting the opposite direction to zero. A site-by-edge matrix was built based on the connection diagram and all edges were weighted using the connectivity values previously calculated with sea current speed and direction. Finally, AEMs were obtained by decomposing the spatial weighted matrix [START_REF] Blanchet | Modelling the effect of directional spatial ecological processes at different scales[END_REF]. All eigenvectors were obtained using the r package 'adespatial' [START_REF] Dray | adespatial: Multivariate multiscale spatial analysis (Version 0.3-8)[END_REF] and were used as predictor variables in independent redundancy analyses (RDA) for each hypothesis. RDA is a constrained ordination analysis, which seeks the linear combination of the predictor variables that best explains the variation of the multivariate response variable. Allelic frequencies for each population were Hellinger-transformed and decomposed by PCA. All the PCs that explain more than 5% of the variation were employed as multivariate response variable in RDAs. From saturated models built for each hypothesis, the eigenvectors that best explained the genetic variation were selected using forward selection with 100,000 permutations. These significant variables were combined into a full model. Partial RDAs and variation partitioning with 100,000 permutations were performed to evaluate the relative contribution of the hypotheses tested. Redundancy analyses and variation partitioning were performed using the 'vegan' r package [START_REF] Oksanen | vegan: Community Ecology package[END_REF].

| RE SULTS

After processing the raw reads with 'ipyrad', 11,624 loci were obtained. Of these, 7,176 were polymorphic, bearing a total of 23,085 SNPs. After removing replicated individuals and poorly sequenced loci, and keeping only biallelic SNPs, there were finally 746 loci and 1,806 SNPs of 240 individuals with 40.46% missing data.

| Genetic diversity, differentiation and population structure

Averaged values across loci and populations of H O , H E , F IS , G ST , θ and Jost's D are shown in Table 1 andTable S1.2. No significant relationship was found between heterozygosity and island size (H O : r = .058, p = .859; H E : r = .231, p = .470). The northern population of Isabela and that of Santa Fe presented the highest value of pairwise differentiation (θ = 0.145 when including the continent; θ = 0.128 within Galápagos), whereas the minimum value was found between the northern population of Isabela and the western population of Santa Cruz (θ = -0.172 when including the continent; θ = -0.109 within Galápagos) (Table S1.3). AMOVA results show a low percentage of variance occurring between islands and between populations within islands, both including (1.46% and 0.91%, respectively) and excluding samples from mainland Ecuador (0.96% and 2.63%, respectively) (Table S1.4).

PCA results showed that none of all the possible combinations of the 11 first principal components (variance explained: 28.38%) grouped the individuals according to geography (Figures S5.1 and S5.2). Neither DAPC nor 'STRUCTURE' results were in agreement with the geographical distribution of the samples and no consensus was found regarding which K-value represents the most plausible number of clusters based on the genetic data (Figure 2a and Figures S5.3-S5.6). DAPC run with groups defined by the user-that is, groups by islands and continent, K = 13-detected highly mixed groups and a clearer genetic identity of individuals from the mainland, Santa Fe and Fernandina (Figure 2b).

| Isolation by resistance, barrier or geographical distance

Resistance by sea currents and low resistance to crossing landmasses (IBR2) was preferred over the competing hypotheses in the reciprocal causal modelling framework based on the relative support values of the matrix, with positive values in the IBR2 column and negative values in the corresponding row (Figure 3a). Furthermore, Mantel tests and some of the partial Mantel tests that include IBR2 and IBD as predictor hypotheses of the genetic distances were statistically significant (Figure 3b). These results did not change when analysing the rainy (from December to May) and dry seasons (from June to November) separately (Tables S1.5 andS1 Redundancy analyses performed with eigenvector maps as predictor variables also support IBR2 as the hypothesis that best explains the spatial genetic variability. Nine PCs retrieved from the PCA of the Hellinger-transformed allele frequencies were used as response variable in the RDAs. None of the dbMEM variables were found to be significant (IBD hypothesis). Two MEM variables (MEM4 and MEM11) significantly explained 9.40% of genetic variability included in the PCs (IBB hypothesis; R 2 a = .094; F = 2.090; p = .040). Two AEMs built under the IBR1 hypothesis (only sea dispersal permitted) were significant (AEM2 and AEM14), explaining 10.54% of the response variability (R 2 a = .105; F = 2.237; p = .003). Four AEMs built under the IBR2 hypothesis (sea and inland dispersal) significantly explained 19.47% of the genetic variance (AEM3, AEM10, AEM11, AEM19; R 2 a = .194; F = 2.269; p = .003) and modelled a spatial genetic structure strongly determined by the prevailing westward sea currents (Figure 4). The full model that combines all significant eigenvectors explained 28.03% of the variance (R 2 a = .2803; F = 2.022; p = .003), but only two AEMs from the IBR2 hypothesis (AEM10 and TA B L E 1 Observed heterozygosity (H O ), expected heterozygosity (H E ) and inbreeding coefficient (F IS ) across loci and populations of Galápagos, mainland Ecuador and populations of Galápagos and the continent considered together. Indices of fixation and differentiation are also shown considering only Galápagos populations and including the continental samples AEM11) and MEM4 (IBB hypothesis) were found to be significant (p < .05). Results of variation partitioning showed that the relative contribution of the AEMs from the IBR2 hypothesis significantly explained 16.86% of the genetic variation (p = .006), while the contribution of IBB eigenvectors was not significant (p = .085) and explained 6.79% of the variance. The overlapping effect of both types of variables was 2.61% and the residual variation 73.74%.

| D ISCUSS I ON

Cryptocarpus pyriformis maintains high gene flow among its populations, which is revealed by low genetic differentiation and little genetic structure, even when including continental populations. This result agrees with the hypothesis of wide dispersal of C. pyriformis through thalassochory, the species being highly admixed at archipelago scale. Sea currents were found to be the main driver of genetic differentiation between islands within the Galápagos, and inland dispersal may have contributed to defining the spatial genetic pattern found within islands.

| Thalassochorous species have high dispersal ability

Plant diaspores adapted to sea dispersal are able to travel across long distances, even thousands of kilometres [START_REF] Cerón-Souza | Contrasting demographic history and gene flow patterns of two mangrove species on either side of the Central American Isthmus[END_REF][START_REF] Miryeganeh | Long-distance dispersal by sea-drifted seeds has maintained the global distribution of Ipomoea pes-caprae subsp. brasiliensis (Convolvulaceae)[END_REF][START_REF] Nathan | Mechanisms of long-distance seed dispersal[END_REF][START_REF] Ridley | The dispersal of plants throughout the world[END_REF][START_REF] Takayama | Gene flow and population subdivision in a pantropical plant with sea-drifted seeds Hibiscus tiliaceus and its allied species: Evidence from microsatellite analyses[END_REF]. This high dispersal ability favours species arrival and colonization of oceanic islands isolated amidst the ocean, confirmed by the fact that thalassochorous species are among the first pioneer plants on newly emerged islands [START_REF] Abe | Colonization of Nishino-Shima island by plants and arthropods 31 years after eruption[END_REF][START_REF] Magnússon | Plant colonization, succession and ecosystem development on Surtsey with reference to neighbouring islands[END_REF][START_REF] Thornton | Krakatau: The destruction and reassembly of an island ecosystem[END_REF]. Thalassochorous species are usually widely distributed across oceanic archipelagos [START_REF] Arjona | Long-distance dispersal syndromes matter: Diaspore-trait effect on shaping plant distribution across the Canary Islands[END_REF][START_REF] Vargas | Plant colonization across the Galápagos Islands: Success of the sea dispersal syndrome[END_REF] and high inter-island gene flow has been reported for some of them, such as Calophyllum inophyllum [START_REF] Hanaoka | Genetic structures of Calophyllum inophyllum L., a tree employing sea-drift seed dispersal in the northern extreme of its distribution[END_REF] and Terminalia catappa [START_REF] Setsuko | Genetic variation of pantropical Terminalia catappa plants with sea-drifted seeds in the Bonin Islands: Suggestions for transplantation guidelines[END_REF] in the Bonin Islands, and Pandanus tectorius in Hawai'i [START_REF] Gallaher | Navigating the 'broad freeway': Ocean currents and inland isolation drive diversification in the Pandanus tectorius complex (Pandanaceae)[END_REF]. In accordance with the high dispersal ability reported for thalassochorous species, the salt bush has been able to colonize Galápagos from the continent (about 1,000 km away), possibly several times independently. The species has spread across all the 13 largest islands of the archipelago and populations maintain high genetic connectivity.

| Oceanic currents as dispersal vector

The high dispersal ability of thalassochorous plants has led many authors to consider them as panmictic species for which no differentiation is expected, considering the ocean as barrier-free for these species [START_REF] Carlquist | The biota of long-distance dispersal. V. Plant dispersal to Pacific Islands[END_REF][START_REF] Nathan | Mechanisms of long-distance seed dispersal[END_REF][START_REF] Paulay | Diversification in the tropical Pacific: Comparisons between marine and terrestrial systems and the importance of founder speciation[END_REF]. However, genetic studies have shown that this statement is not ubiquitous, as some level of genetic structuring has been found in sea-dispersed plants, not only mediated by land barriers but because of sea current properties, in particular current speed and direction [START_REF] Cerón-Souza | Contrasting demographic history and gene flow patterns of two mangrove species on either side of the Central American Isthmus[END_REF][START_REF] Gallaher | Navigating the 'broad freeway': Ocean currents and inland isolation drive diversification in the Pandanus tectorius complex (Pandanaceae)[END_REF][START_REF] Guo | Genetic discontinuities in a dominant mangrove Rhizophora apiculata (Rhizophoraceae) in the Indo-Malesian region[END_REF][START_REF] Hanaoka | Genetic structures of Calophyllum inophyllum L., a tree employing sea-drift seed dispersal in the northern extreme of its distribution[END_REF]; Hodel [START_REF] Ngeve | Contrasting effects of historical sea level rise and contemporary ocean currents on regional gene flow of Rhizophora racemosa in eastern Atlantic mangroves[END_REF][START_REF] Wee | Oceanic currents, not land masses, maintain the genetic structure of the mangrove Rhizophora mucronata Lam. (Rhizophoraceae) in Southeast Asia[END_REF]. Despite the importance of sea currents in long-distance dispersal, few studies have statistically assessed their role in genetic differentiation and distribution of terrestrial organisms. Nevertheless, several authors [START_REF] Fajardo | The role of abiotic mechanisms of long-distance dispersal in the American origin of the Galápagos flora[END_REF][START_REF] Harvey | Spatial patterns of inter-island plant and bird species movements in the Galapagos Islands[END_REF] have evaluated the relationship between sea currents and species distributions, and also, for instance, the relationship between sea currents and gene flow of Caribbean lizards [START_REF] Calsbeek | Ocean currents mediate evolution in island lizards[END_REF] and mangroves [START_REF] Hodel | Terrestrial species adapted to sea dispersal: Differences in propagule dispersal of two Caribbean mangroves[END_REF]. All of these found a significant contribution of sea currents to the distribution and genetic patterns. In agreement with these studies, genetic differentiation of the Galápagos salt bush populations is correlated to direction and speed of sea currents, as shown by reciprocal causal modelling and redundancy analysis results. This is as expected for a sea-dispersed species.

In addition, geographical and genetic distances were revealed to be significantly correlated by Mantel test (Figure 3b), as cost-distances increase not only with the resistance to dispersal driven by oceanic currents, but geographical distance as well. Similar results were found in mangroves by [START_REF] Hodel | Terrestrial species adapted to sea dispersal: Differences in propagule dispersal of two Caribbean mangroves[END_REF].

In addition to sea dispersal, our results indicate that inland gene flow is likely at work for this coastal species, as only the hypothesis that combines resistance by maritime and terrestrial dispersal (IBR2)

was consistently selected by both reciprocal causal modelling and redundancy analysis. Finches and lizards that consume salt bush fruits might act as potential inland dispersers of this species [START_REF] Heleno | Frugivory and seed dispersal in the Galápagos: What is the state of the art?[END_REF][START_REF] Heleno | Seed dispersal networks in the Galápagos and the consequences of alien plant invasions[END_REF]. However, analysis of genetic data mainly from the nuclear genome prevents distinguishing between the seed and pollen flow contributions to such short-distance gene flow. Indeed, insect-mediated pollen dispersal has been proposed to explain the connectivity found across the Panama Isthmus between populations of Hibiscus tiliaceus [START_REF] Takayama | Gene flow and population subdivision in a pantropical plant with sea-drifted seeds Hibiscus tiliaceus and its allied species: Evidence from microsatellite analyses[END_REF] and Ipomoea pescaprae [START_REF] Miryeganeh | Long-distance dispersal by sea-drifted seeds has maintained the global distribution of Ipomoea pes-caprae subsp. brasiliensis (Convolvulaceae)[END_REF].

| Quaternary climatic and sea-level changes affected coastal plant populations

Historical geological changes-such as the emergence and subsidence of oceanic islands (e.g. [START_REF] Christie | Drowned islands downstream from the Galapagos hotspot imply extended speciation times[END_REF])-together with climatic and sea-level changes might have altered oceanic currents [START_REF] Gillespie | Long-distance dispersal: A framework for hypothesis testing[END_REF]. In addition, coastal habitats are expected to have suffered frequent disturbance due to past climate changes and variation in sea level, causing successive episodes of flooding, drought and other extreme meteorological events, possibly leading to local extinctions [START_REF] Kudoh | Consequences of hydrochory in Hibiscus[END_REF][START_REF] Ngeve | Contrasting effects of historical sea level rise and contemporary ocean currents on regional gene flow of Rhizophora racemosa in eastern Atlantic mangroves[END_REF] and further recolonizations [START_REF] Ali | Exploring the combined role of eustasy and oceanic island thermal subsidence in shaping biodiversity on the Galápagos[END_REF]. Such metapopulation dynamics have probably been dramatically catalysed by active volcanism, which has transformed the landscape of many of the Galápagos Islands for more than two million years since their emergence [START_REF] Geist | Paleogeography of the Galápagos Islands and biogeographical implications[END_REF]. Volcanism and sea-level changes may have had an impact on the demographic history of coastal plant species, likely causing profound changes in their effective population sizes through time. As a consequence, historical processes are expected to have left an imprint on the gene flow detected among populations of the Galápagos salt bush, almost ruling out the possibility of finding any stronger relationship with the current ocean circulation [START_REF] Hodel | Terrestrial species adapted to sea dispersal: Differences in propagule dispersal of two Caribbean mangroves[END_REF]. In addition to the weak genetic structure found for the Galápagos populations in this study, the low heterozygosity and its

  tween populations: (a) isolation by resistance due to ocean currents with no overland dispersal possible (IBR1); (b) isolation by resistance due to ocean currents and uniform low resistance in crossing landmasses (IBR2); (c) isolation by barrier (IBB) with islands being isolated due to the impossibility of crossing the sea barrier; and (d) isolation by geographical distance (IBD). The net contribution of the hypothesis components to the overall genetic variation was assessed using Mantel and partial Mantel tests in a reciprocal causal modelling framework and, independently, by decomposing the different distance matrices into eigenvector maps and using them as predictors in redundancy analyses.

  .6). Besides this, no significant correlation was found between monthly El Niño Southern Oscillation index and monthly Mantel's r (r M ) comparing genetic and cost distances (IBR2), indicating a negligible effect of ocean current oscillations due to El Niño events (Figure S4.1). The isolation-by-distance (IBD) and isolation-by-resistance hypotheses (IBR1 and IBR2) were highly correlated (IBR1 ~ IBR2: r M = .498, p < .001; IBD ~ IBR1: r M = .787, p < .001; IBD ~ IBR2: r M = .757; p < .001).

F I G U R E 2

 2 Results of the discriminant analysis of principal components (DAPC) for the genetic structure of Cryptocarpus pyriformis across the sampled area (Galápagos and mainland Ecuador). Compoplots (bar graphs) show the genetic composition of each population and scatter plots show how similar the established groups are. (a) K = 5 genetically similar groups established using the k-means algorithm. (b) Clusters determined using populations a priori (K = 13)

F

  I G U R E 3 Reciprocal causal modelling results. (a) Hypothesis testing. The coefficient in each cell is the result of the difference between Mantel's r (r M ) of two partial Mantel tests: (i) effect of the focal hypothesis (column) partialling out the effect of an alternative hypothesis (row), and (ii) effect of the alternative hypothesis partialling out the effect of the focal hypothesis (focal|alternative -alternative|focal). The colour gradient from blue to red varies according to the coefficient in the cell and indicates the relative support of the focal hypothesis. The best supported hypothesis is IBR2, as all the coefficients in the IBR2 column are positive, and negative in the corresponding row. (b) Mantel's r for simple Mantel tests (diagonal) and partial Mantel tests that evaluate the effect on the genetic differentiation of a focal hypothesis (column) controlling for the effects of an alternative hypothesis (row). Statistically significant correlations (p < .05) are shown by * 2018;[START_REF] Miryeganeh | Long-distance dispersal by sea-drifted seeds has maintained the global distribution of Ipomoea pes-caprae subsp. brasiliensis (Convolvulaceae)[END_REF][START_REF] Ngeve | Contrasting effects of historical sea level rise and contemporary ocean currents on regional gene flow of Rhizophora racemosa in eastern Atlantic mangroves[END_REF]

  Cryptocarpus pyriformis is abundant in coastal xeric areas of the archipelago. This pattern favours short inland pollen dispersal by insects and birds, both effective pollinators of the species (Hervías-Parejo & Traveset, 2018; Olesen et al., 2018). The respective contribution of pollen versus seed dispersal in explaining within-islands isolation by distance remains an open question.

  lack of relationship with island size could indicate a recent colonization or recolonization of the archipelago, or at least that the populations underwent similar demographic processes in the recent past. Further avenues for research are opened to investigate the demographic history of the Galápagos salt bush and the clues it provides regarding the geological history of the archipelago. Although some open questions remain, our study highlights the usefulness of combining methods developed under a joint landscape and seascape genetic framework, to statistically assess how the F I G U R E 4 Redundancy analysis (RDA) biplot showing the two first canonical axes (RDA1 and RDA2 -percentage of total variation explained shown in brackets) of the canonical ordination of the four significant asymmetric eigenvector maps (AEM variables-arrows) built under the hypothesis that considers together maritime and terrestrial dispersal (IBR2), and the 23 Galápagos populations (dots coloured according to geographical position-see map-legend). Island abbreviations: GE = Genovesa, MA = Marchena, PT = Pinta, SF = Santa Fe, SX = Santa Cruz, SA = Santiago, PZ = Pinzón, IS = Isabela, FE = Fernandina, SC = San Cristóbal, ES = Española, FL = Floreana surrounding environment modulates gene flow between populations of coastal plants.

ACK N OWLED G EM ENTS

We thank Emilio Cano and Morgan Gueuning for their help in the wet lab and all N.A.'s group for their advice on technical issues, Raphaël Leblois for his advice on population genetic analyses, Javier Morente for useful discussion about the manuscript, and two anonymous referees, the Chief and Handling Editors for their valuable comments that helped to improve the manuscript.

We also thank the Galápagos fieldwork team, the Charles Darwin Foundation and the Galápagos National Park, which facilitated plant material sampling (permit numbers: PC-46-14 and PC-21-15). YA was funded by the Ministerio de Economía y Competitividad (Spain) (BES-2013-062843; EEBB-I-15-10273; EEBB-I-16-11437). This study is framed within a project also financed by the Ministerio de Economía y Competitividad (CGL2012-C02-01). The manuscript was edited by Guido Jones, currently funded by the Cabildo de Tenerife, under the TFinnova Programme supported by MEDI and FDCAN funds.

/Award Number: BES-2013-062843, CGL2012-C02-01, EEBB-I-15-10273 and EEBB-I-

DATA AVA I L A B I L I T Y S TAT E M E N T

Final assembled and filtered SNP dataset and scripts are stored and publicly accessible in figshare repository (https://doi.org/10.6084/ m9.figsh are.12529 523.v1). Modified rWind functions to work with sea current data can be found in GitHub repository (https:// github.com/jabio logo/rWind/ tree/maste r/man). Raw demultiplexed sequences of the ddRADseq dataset are available on NCBI SRA (BioProject Accession PRJNA655805).

S U PP O RTI N G I N FO R M ATI O N

Additional supporting information may be found online in the Supporting Information section.

How to cite this article: Arjona Y, Fernández-López J, Navascués M, Alvarez N, Nogales M, Vargas P. Linking seascape with landscape genetics: Oceanic currents favour colonization across the Galápagos Islands by a coastal plant. J Biogeogr. 2020;00:1-12. https://doi.org/10.1111/jbi.13967