Christophe Vieira

Labex-Corail Coréus

Sofie D'hondt

Olivier De Clerck

Claude E Payri

TOWARD AN INORDINATE FONDNESS FOR STARS, BEETLES AND LOBOPHORA? SPECIES DIVERSITY OF THE GENUS LOBOPHORA (DICTYOTALES, PHAEOPHYCEAE) IN NEW CALEDONIA 1

Keywords: Lobophora, GMYC, ABGD, species delimitation, New Caledonia, new species, phylogeny, taxonomy 4 ABGD, automated barcoding gap discovery, AIC, Akaike information criterion, ANOVA, analysis of variance, BI, Bayesian inference, bGMYC, Bayesian implementation of the general mixed Yule coalescent, CTAB, cetyltrimethyl ammonium bromide, GMYC: general mixed Yule coalescent, GTR, generalized time reversible, MCMC molecular operational taxonomic unit, Tukey

Until the recent use of molecular markers, species diversity of Lobophora, an ecologically important brown algal genus with a worldwide distribution in temperate and tropical seas, has been critically underestimated. Using a DNA-based taxonomic approach, we re-examined diversity of the genus from New Caledonia in the Southwest Pacific Ocean. First, species were delineated using GMYC-based and barcoding gap approaches applied to a mitochondrial cox3 dataset. Results were subsequently confirmed using chloroplast psbA and rbcL datasets. Species delimitation analyses agreed well across markers and delimitation algorithms, with the barcoding gap approach being slightly more conservative. Analyses of the cox3 dataset resulted in 31 to 39 molecular operational taxonomic units, four of which are previously described species (L. asiatica , L. crassa, L. nigrescens s.l., L. pachyventera). Of the remaining MOTUs for which we obtained a representative number of sequences and results are corroborated across analyses and genes, we describe ten species de novo: L. abaculusa, L. abscondita, L. densa, L. dimorpha, L. gibbera, L. hederacea, L. monticola, L. petila, L. rosacea, and L. undulata. Our study presents a excellent case of how a traditional morphology-based taxonomy fails to provide accurate estimates of algal diversity. Furthermore, the level of Lobophora diversity unveiled from a single locality in the Pacific Ocean raises important questions with respect to the global diversity of the genus, the distributions and range sizes of the individual species, as well as the mechanisms facilitating co-existence.

Introduction

Contrary to substantial historical disagreement on the generic classification of the genus Lobophora J.Agardh (J.V.Lamouroux 1809, C.Agardh 1817, J. [START_REF] Agardh | Analecta algologica, observationes de speciebus algarum minus cognitae earumque dispositione: Continuatio I. . Lunds Universitets Års-Skrift, Andra[END_REF][START_REF] Papenfuss | Notes on algal nomenclature[END_REF][START_REF] Womersley | A critical survey of the marine algae of southern Australia. II. Phaeophyta[END_REF]), species-level taxonomy has been remarkably stable. Traditionally only three Lobophora species were recognized, with L. variegata (J.V.Lamouroux) Womersley ex E.C.Oliveira being by far the most commonly reported species. Literature data make it seem that L. variegata is widely distributed in temperate to tropical parts of the Atlantic (incl.

Mediterranean Sea), Indian and Pacific Ocean. The other two species L. papenfussii (W.R.Taylor) Farghaly and L. dichotoma (R.H.Simons) P.C.Silva were only sporadically reported from the Indo-Pacific and South Africa respectively. From 2000 until 2012, three more species were described (L. minima V. Krishnamurthy and M.Baluswami (2000), L. indica V. Krishnamurthy and M.Baluswami (2000) and L. rickeri [START_REF] Kraft | Algae of Australia: Marine Benthic Algae of Lord Howe Island and the Southern Great Barrier Reef, 2: Brown algae[END_REF]), based on morphological criteria only.

From a molecular phylogenetic perspective Lobophora had not received much attention (but see [START_REF] Hoshina | Molecular phylogeny of the Dictyotaceae (Phaeophyceae) with emphasis on their morphology and its taxonomic implication[END_REF][START_REF] Phillips | Resolving evolutionary relationships among the brown algae using chloroplast and nuclear genes[END_REF][START_REF] Bittner | Molecular phylogeny of the Dictyotales and their position within the Phaeophyceae, based on nuclear, plastid and mitochondrial DNA sequence data[END_REF]) until a recent study of [START_REF] Sun | Taxonomic revision of the genus Lobophora (Dictyotales, Phaeophyceae) based on morphological evidence and analyses rbc L and cox3 gene sequences[END_REF]. The latter authors recognized nine major Lobophora clades based on chloroplast rbcL and mitochondrial cox3 gene sequences, four of which were formally described as new species (i.e. L. asiatica Z.Sun, Ji.Tanaka and H.Kawai, L. crassa Z.Sun, P.-E.Lim and H.Kawai, L. pachyventera Z.Sun, P.-E.Lim, Tanaka and H.Kawai, L. australis Z. Sun, Gurgel and H.Kawai).

In total, 10 species are currently accepted taxonomically (Guiry and Guiry, 2013).

Despite the ecological importance of Lobophora in seaweed-coral-grazing interactions and competition (De Ruyter van Steveninck and Breeman 1987a,b, De Ruyter van Steveninck et al. 1988a,b,c, Coen and Tanner 1989[START_REF] Diaz-Pulido | Doom and boom on a resilient reef: climate change, algal overgrowth and coral recovery[END_REF][START_REF] Rasher | Seaweed allelopathy degrades the resilience and function of coral reefs[END_REF][START_REF] Anthony | Ocean acidification and warming will lower coral reef resilience[END_REF][START_REF] Slattery | Allelopathy in the tropical alga Lobophora variegata (Phaeophyceae): mechanistic basis for a phase shift on mesophotic coral reefs?[END_REF], the species diversity of the genus remains largely unaddressed.

Here we study the diversity of Lobophora in New Caledonia. New Caledonia is located just south of the coral triangle, recognized as the global center of marine biodiversity, and displays tropical to subtropical-temperate conditions. The Lobophora flora has been comprehensively sampled over the last decades from various regions and the large amount of material revealed a large morphological diversity associated to the ecological variation justifying the present study.

The paper of [START_REF] Sun | Taxonomic revision of the genus Lobophora (Dictyotales, Phaeophyceae) based on morphological evidence and analyses rbc L and cox3 gene sequences[END_REF] provided two important insights about the genus Lobophora, (1) the existence of a rich and yet to be discovered diversity and (2) the occurrence of cryptic diversity lacking distinctive morphological features between taxa.

Decisions on species concepts as well as the practical criteria to delimit species represent critical aspects for studies aiming to elucidate species level diversity (e.g. [START_REF] Harrison | Linking evolutionary pattern and process[END_REF][START_REF] Agapow | The impact of species concept on biodiversity studies[END_REF]). For algae it has long been recognized that diversity is often inadequately reflected in the organism's morphology. It is therefore not surprising that, coinciding with a growing ease to obtain molecular data, the latter have become the standard for delimiting algal species (see Alverson 2008;[START_REF] De Clerck | Algal taxonomy: a road to nowhere?[END_REF][START_REF] Leliaert | DNA-based species delimitation in algae[END_REF]. Accompanying a growing dependency on DNA sequence data in biodiversity assessment, a variety of approaches and algorithms have been proposed to detect discontinuities in genetic variation representative for species boundaries (e.g. [START_REF] Wiens | Delimiting species using DNA and morphological variation and discordant species limits in spiny lizards (Sceloporus)[END_REF][START_REF] Sites | Operational criteria for delimiting species[END_REF][START_REF] Carstens | How to fail at species delimitation[END_REF].

Since, species delimitation may be influenced by the gene information content as well as the species delimitation method, we test species boundaries in Lobophora using three species delimitation methods, a General Mixed Yule Coalescent (GMYC) model (Pons et al. 2006, Fujisawa and[START_REF] Fujisawa | Delimiting Species Using Single-locus Data and the Generalized Mixed Yule Coalescent (GMYC) Approach: A Revised Method and Evaluation on Simulated Datasets[END_REF], the Bayesian implementation of the GMYC model [START_REF] Reid | Phylogenetic estimation error can decrease the accuracy of species delimitation: a Bayesian implementation of the general mixed Yule-coalescent model[END_REF]) and an Automated Barcoding Gap Discovery method (ABGD) (Puillandre et al. 2011). The combination of several molecular methods for species delimitation is becoming a reference to detect species boundaries and have been used in different taxonomical groups (Jörger et al. 2012 for sea slugs; Kekkoken and Hebert 2014 for moths; [START_REF] Cornils | Evidence of cryptic and pseudocryptic speciation in the Paracalanus parvus species complex (Crustacea, Copepoda, Calanoida)[END_REF] for copepods; [START_REF] Alò | Diversity of Aplochiton Fishes (Galaxiidea) and the Taxonomic Resurrection of A. marinus[END_REF] for fishes). To our knowledge it is the first time that such a combination is used for algae species delimitation.

Species delimitation is in the first place carried out using a mitochondrial cox3 dataset for which we had the most complete taxon sampling. To investigate up to which extent results were influenced by marker choice, analyses were repeated for chloroplast rbcL and psbA datasets, which contained less sequences per taxon compared to the cox3 dataset. Subsequently, we studied the morphology and ecology of the New Caledonian specimens to determine up to which extent the DNA-based species are morphologically and ecologically diverged. Sampling was carried out mainly by SCUBA from 3 down to 90m deep or by snorkeling and reef walking. The specimens were readily stored in a cooler and desiccated in silica gel for subsequent DNA extraction once at the laboratory. Specimens were dried and mounted on herbarium sheets and deposited at the IRD Herbarium of Nouméa (New Caledonia, IRD-NOU).

Materials and Methods

Sampling

For the earliest collections, dry Herbarium specimens were used as DNA source. The New Caledonia samples were complemented with a few collections from Papua New Guinea (Madang 2012) and the Maldive Islands (2011). The origin of the specimens and accession numbers are detailed in TABLE S1.

DNA extraction, amplification, sequencing and phylogenetic analyses

Total genomic DNA was extracted from 235 Lobophora samples, 228 from New Caledonia, 5 from Papua New Guinea and 2 from the Maldive Islands using a CTAB-extraction method [START_REF] De Clerck | A revised classification of the Dictyoteae (Dictyotales, Phaeophyceae) based on rbcl and 26s ribosomal DNA sequence analyses1[END_REF]. Genomic DNA was subsequently purified with a Wizard® DNA Clean-Up System (Promega Inc., Madison, WI, USA) following the manufacturer's instructions.

Sequences were generated from one mitochondrial gene (cox3), two chloroplast genes (psbA, rbcL) and the 5'-end of the nuclear encoded large subunit rDNA (LSU, ca. 1200 bp). PCR and sequencing conditions are detailed in TABLE S2. LSU sequences were not tested for species delimitation because of the low number of sequences obtained, but were integrated in the concatenated alignment to generate a species tree with improved resolution. In addition to the sequences generated in the present study, 25 cox3, 4 psbA, 33 rbcL and 6 LSU Lobophora sequences from GenBank were added to the alignments (TABLE S1). Sequences were aligned using MUSCLE implemented in eBioX 1.5.1 (www.ebioinformatics.org). Ambiguously aligned regions in the LSU alignment were removed by eye.

Species delimitation

Following exploratory ML and Bayesian analyses (results available upon request), ultrametric gene trees were constructed using Bayesian analyses in BEAST v1.7.5 [START_REF] Drummond | Bayesian phylogenetics with BEAUti and the BEAST 1.7[END_REF] for the cox3, rbcL and psbA alignments. A GTR+G substitution model was identified as the bestfitting model for each individual gene, based on the Akaike Information Criterion (AIC) using jModelTest 2 [START_REF] Darriba | jModelTest 2: more models, new heuristics and parallel computing[END_REF]. BEAST analyses were run under a strict molecular clock in combination with a Constant Coalescent tree prior. Other priors were set to default. In order to check for convergence of the MCMC chains, we performed two independent runs for 10 7 generations each, starting from random trees and sampling every 10 4 generations. MCMC output files of the independent runs were inspected in Tracer 1.5 (Rambaut and Drummond 2009) for acceptable effective sample sizes (ESS > 200). A burn-in was applied once log-likelihood values had stabilized. Maximum clade credibility trees and posterior probability for the nodes were calculated using the postburnin trees using TreeAnnotator 1.6.2 (included in the BEAST package). All tree searches were conducted on the Cipres web portal [START_REF] Miller | Creating the CIPRES Science Gateway for inference of large phylogenetic trees[END_REF].

We used a Maximum Likelihood (GMYC) as well as a Bayesian Implementation (bGMYC) of the GMYC model [START_REF] Pons | Sequence-based species delimitation for the DNA taxonomy of undescribed insects[END_REF][START_REF] Reid | Phylogenetic estimation error can decrease the accuracy of species delimitation: a Bayesian implementation of the general mixed Yule-coalescent model[END_REF]. Both methods are able to discriminate between population and speciation patterns on a given ultrametric tree. GMYC analyses under a single-threshold were conducted in R (R Core Team, 2014) using the package "Splits". The bGMYC model was performed using "bGMYC" [START_REF] Reid | Phylogenetic estimation error can decrease the accuracy of species delimitation: a Bayesian implementation of the general mixed Yule-coalescent model[END_REF]) in R using a subsample of 100 trees from the posterior distribution of BEAST as suggested by the authors. Markov chain Monte Carlo (MCMC) chains were run for each tree for 10,000 generations with a burn-in comprising the first 1,000 generations once the log-likelihood values had stabilized, and sampling every 100 generations.

Automatic Barcode Gap Discovery (ABGD, [START_REF] Puillandre | ABGD, Automatic Barcode Gap Discovery for primary species delimitation[END_REF]) is an exploratory tool based on pairwise distances to detect automatically significant difference in intra and inter specific variation (i.e. barcoding gap), without an a priori species hypothesis. These analyses were performed on the abgd website (wwwabi.snv.jussieu.fr/public/abgd/abgdweb.html. Accessed 2013 October 12) selecting default parameters except for the relative gap width (X) which was set to 1 and the number of steps which was set to 100. The distance matrix was build under a K2P model. Species boundaries were subsequently defined based on the congruence of the three methods and are detailed in the discussion.

Species tree inference

Maximum likelihood (ML) and Bayesian inference (BI) phylogenetic species trees were generated from a concatenated alignment including cox3 (610 bp), psbA (919 bp), rbcL (1360 bp) and LSU rDNA (1361 bp) genes, partitioned by gene and codon position. The concatenated alignment contained a single representative per Molecular Operational Taxonomic Unit (MOTU) resulting from the species delineation analyses of the rbcL dataset. The matrix was 70% filled at the MOTU level. A selection of Zonaria C.Agardh (Dictyotales, Phaeophycea), Padina Adanson (Dictyotales, Phaeophycea) and Dictyota J.V.Lamouroux (Dictyotales, Phaeophycea) species were used as outgroup taxa (cf. TABLE S1). ML analyses were conducted using RAxML under a GTR+CAT model [START_REF] Stamatakis | RAxML-VI-HPC: Maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models[END_REF]. The robustness of the resulting phylogenies was tested using 1000 replicates of a rapid bootstrap heuristic [START_REF] Stamatakis | A rapid bootstrap algorithm for the RAxML web servers[END_REF]. BI, using MrBayes v3.2.2 [START_REF] Ronquist | MrBayes 3: Bayesian phylogenetic inference under mixed models[END_REF], initiated with a random starting tree and ran four chains of MCMC iterations simultaneously for 100 million generations. The first 100,000 (25%) trees sampled were discarded as burn-in, based on the stationarity of lnL as assessed using Tracer version 1.5 (Rambaut and Drummond 2009). A consensus topology and posterior probability values were calculated from the remaining trees.

Morphological and ecological analyses

Morphological observations of Lobophora species included analyses of the external and internal (anatomy) structure of the specimens. Based on our field observations we distinguished the occurrence of seven main growth forms, namely (1) stipitate, (2) fasciculate, (3) conk-like, (4) decumbent, (5) anastomosing, (6) procumbent and (7) crustose as illustrated and defined in FIG.

2. For the internal morphology, longitudinal and transverse sections were made of the basal, middle and distal portions of the thallus using a medical freezing portable microtome (Labonord®). Photographs of the sections were taken with a digital camera (Olympus Camedia C-5050 5.0 Megapixel, Tokyo, Japan) attached to a compound microscope (Olympus BH-2, Tokyo, Japan). The number and size of the cortical (dorsal and ventral) and medulla cells of the basal, middle and distal portions of the thallus were measured as shown in FIG. 3, which resulted in the measurements of 9 anatomical traits (i.e. number of dorsal and ventral cells; total number of cells; thallus thickness; dorsal, medullar and ventral heights; medullar width and length). The surface of the thallus with rhizoids was defined as the ventral surface. A total of 285 specimens, from one to 15 specimens per species, were examined for morphological analyses. Every specimen studied morphologically has been sequenced for at least the cox3 marker. A few sequences which were too short were not included in the molecular analyses. Descriptive statistics were generated for the anatomical traits and correlations between them were tested to select independent traits for subsequent univariate analyses. Mean anatomical traits were tested for equality by a one-way ANOVA and post-ANOVA Tukey Honestly Significant Difference (HSD) tests. The data were tested for normality and homogeneity of variances by means of a Shapiro-Wilk test of normality and the Bartlett test of the homogeneity of variances. The thickness data were log-transformed prior to analysis, to meet assumptions of normality and homogeneity of variance. All analyses were conducted using R. Ecologically, we identified three major substratum preferences in the field specific to some groups of species: (1) niched among or growing on live corals, (2) growing at the base of live corals, on dead corals, coral rubbles or bedrock and (3) growing niched among Sargassum beds.

Results

Species delimitation

Species delimitation based on the cox3 alignment (610 bp x 210 sequences) using GMYC under a single threshold resulted in an estimate of 37 for MOTUs, with a confidence interval of 36-49 (FIG. 4). The number of specimens per MOTU ranged from 1 (singletons) to 45 with an average of 6.5. bGMYC analysis of posterior probabilities of conspecificity within cox3 Lobophora clusters was high (P > 0.9) and resulted in a species delimitation which was marginally less conservative than GMYC, differing in 2 cases only (FIG. 4): IRD10187 was resolved as a singleton (prob. 0.59), d271 and d6625 were resolved as a separate cluster (prob. 0.648). The ABGD approach is slightly more conservative, grouping four MOTUs that were split in both GMYC analyses (FIG. 4). Species delimitation analyses were repeated for rbcL (1345 bp x 139 sequences) and psbA (919 bp x 88 sequences) datasets to investigate if the cox3 results were stable across genes. In all analyses the likelihood of the GMYC model was significantly higher (p < 0.001) than that of the null model of uniform coalescent branching rates. GMYC analyses of rbcL data yielded (40-) 47 13 (-54) MOTUs while the psbA data resulted in (17-) 19 (-34) MOTUs (TABLE 1). Contrary to the cox3 dataset, no incongruence between the various delimitation methods was detected.

Unequal sampling across markers complicates a detailed comparison of results from different markers, but even without a fully congruent sampling it was clear that the outcome of the analyses was stable across genes (TABLE S3). Six MOTUs from the cox3 ABGD analysis were subdivided in less inconclusive units in the rbcL dataset. All but one of the cox3 bGMYC MOTUs on the other hand were confirmed in the rbcL dataset. Data from the psbA dataset are less informative because of the high number of missing MOTUs (47%), but of the cox3 ABGD MOTUs present 2 are subdivided and one is merged with another MOTU. Similarly, 2 cox3 bGMYC MOTUs are merged. In addition, inclusion of Genbank accessions in the rbcL dataset yielded 9 additional MOTUs, which were not represented in either the cox3 or psbA dataset. This resulted in the rbcL gene alignment being the most diverse in terms of MOTUs, but with a significantly higher number of singletons than cox3.

Morphological and ecological characters

The morphology and ecology of the specimens from New Caledonia were studied to determine up to which extent the MOTUs are morphologically and ecologically diverged. For practical reasons we introduce names of newly described species already in the sections below. Results and interpretations of correlation analyses between the nine anatomical characters measured are given in the supplementary text (TABLE S4). Boxplots were used to show inter-and intraspecific variation of six anatomical traits (thallus thickness; dorsal and ventral height; medulla height, width and length) (FIG. 5). Anatomical characters related to cell height differed significantly among species as well. On the other hand cell length and width displayed some variation but were overall less diagnostic. Among the three independent anatomical traits (i.e. thallus thickness, medulla width and length), the thallus thickness presented the most significant interspecific variability and was therefore retained as the only variable for the ANOVA analysis.

The thallus thickness ranged from an average of 57 µm for the thinnest species (L. petila) to 407 µm for the thickest species (L. densa). A continuous grade from these two extreme values was observed and the thickness of several species overlapped. The amount of intraspecific variation differed, with the thicker species presenting a greater variability. A one-way ANOVA analysis (TABLE S4) revealed statistically significant differences and subsequent post-hoc analyses (Tuckey HSD) (TABLE S5) confirmed significant difference between the species thallus thickness means. Seven species presented unique means and distribution (L. densa, L. crassa, L. gibbera, L. hederacea, L. monticola, L. pachyventera and L. petila) and three groups of species exhibited neighboring mean values with comparable variances (FIG. S1). Consequently, thallus thickness may serve to identify seven New Caledonian species but for some groups of species does not suffice to go down to the species level delineated with the phylogenetic approaches.

However, for those 3 groups with similar thickness, external morphology and ecology allow species differentiation (see below).

Species phylogeny

ML and BI analyses of the concatenated alignment (cox3 + rbcL + psbA + LSU) including every MOTU discovered in the species delimitation analyses, yielded similar tree topologies except for the relationships between the MOTUs 29 to 32, and the MOTUs 45 to 47. Results are presented using the BEAST ultrametric tree topology (FIG. 6). The 4-genes analyses resulted in a fairly well-resolved phylogeny with moderate to strong support for most nodes. The phylogenetic tree revealed 6 well-supported lineages (defined as a sequence of species or MOTUs; Lineage A-F) (FIG. 6). However, the position of the MOTU 46 from Guadeloupe, for which only the rbcL sequence is available, is incongruent between the trees. In the BEAST and ML trees MOTU 46 is part of the Lineage A (FIG. 6 andS2), while it comes outside of the Lineage A, in the most basal position, in the Bayesian tree (FIG. S3). This inconsistency may be resolved by acquiring extra sequences for the missing markers, and for the time being we will consider it as part of the Lineage A.

Discussion

Species diversity and taxonomy

In this study we aimed to characterize the diversity of the genus Lobophora in New Caledonia in the South West Pacific Ocean and subsequently address the evolutionary relationships of the New Caledonian representatives. Thereto, we applied the most comprehensive sampling of the genus to date. Although expecting some levels of cryptic or pseudocryptic diversity, much to our initial astonishment cox3 species delimitation analyses yielded between 31, 37 and 39 MOTUs based on ABGD, GMYC and bGMYC analyses, respectively. Both GMYC-based methods were highly congruent. The bGMYC analyses segregated one specimen (IRD10187) from Lobophora crassa2. Likewise, d271 and d6625 were segregated from L. nigrescens s.l. Both results, however, were only moderately supported in the bGMYC analysis, with posterior probabilities of 0.591 and 0.648 respectively. The barcoding gap method yielded a more conservative estimate, but most discrepancies were limited to the L. crassa and L. pachyventera complexes as defined by [START_REF] Sun | Taxonomic revision of the genus Lobophora (Dictyotales, Phaeophyceae) based on morphological evidence and analyses rbc L and cox3 gene sequences[END_REF] and discussed below.

Subsequent analyses of rbcL and psbA dataset were highly congruent with the GMYC and bGMYC results and indicated that the ABGD estimate of the cox3 dataset is likely somewhat over-conservative (TABLE 2 and TABLE S6). Possibly the small sample size of some MOTUs may result in larger units as identified by the barcoding gap approach [START_REF] Jorger | Barcoding against a paradox? Combined molecular species delineations reveal multiple cryptic lineages in elusive meiofaunal sea slugs[END_REF][START_REF] Puillandre | ABGD, Automatic Barcode Gap Discovery for primary species delimitation[END_REF]). We identified one case in which the cox3 GMYC analyses were too conservative (SAP109520) compared to rbcL results, and one case in which they were too liberal (IRD10187). In both situation a single specimens was either added to or segregated from a MOTU.

Our analyses disclosed the occurrence of 29 MOTUs in New Caledonian. These results confirm findings by [START_REF] Sun | Taxonomic revision of the genus Lobophora (Dictyotales, Phaeophyceae) based on morphological evidence and analyses rbc L and cox3 gene sequences[END_REF] of undescribed species diversity in Lobophora. Species boundaries as defined by [START_REF] Sun | Taxonomic revision of the genus Lobophora (Dictyotales, Phaeophyceae) based on morphological evidence and analyses rbc L and cox3 gene sequences[END_REF] of L. asiatica, L. nigrescens sensu [START_REF] Sun | Taxonomic revision of the genus Lobophora (Dictyotales, Phaeophyceae) based on morphological evidence and analyses rbc L and cox3 gene sequences[END_REF] (subsequently referred to as L. nigrescens s.l.) and L. australis are mirrored by our species delimitation.

However, their species delineation appeared to be more conservative for L. crassa and L. pachyventera. GMYC and bGMYC analyses split the L. crassa and L. pachyventera complexes into five and four MOTUs respectively for cox3 (FIG. 4). In the L. crassa complex the New Caledonian specimens were resolved as separate MOTUs, L. crassa2, L. crassa4 and L. crassa5.

Likewise, in the L. pachyventera complex the New Caledonian specimens were resolved as a separate MOTU, L. pachyventera2. However, it should be noticed that the cox3 ABGD results group the L. crassa MOTUs and L. pachyventera1, L. pachyventera2 and L. pachyventera3 in two clusters only. Four of the New Caledonian MOTUs, were assigned to existing species or species complexes (Lobophora crassa, Lobophora asiatica, Lobophora pachyventera and L. nigrescens s.l.). In addition, none of our samples matched the descriptions of the four Lobophora species for which no molecular data are available (i.e. L. variegata, L. dichotoma, L. rickeri, L. papenfussii). The remaining MOTUs could therefore qualify as putative species.

Decisions as to which of these putative new species should be described de novo are based on the availability of a representative set of specimens for a single MOTU and congruence between the various species delimitation algorithms. In this we opt for a conservative approach, describing only those species for which we had (1) at least 3 sequences (specimens) for cox3, (2) at least sequences for the three markers (cox3, rbcL and psbA), and (3) which resulted in consensual results between analyses (GMYC, bGMYC and ABGD) and genes. In other words, we opted for the least inclusive species delimitation. Based on this rationale we describe 10 species de novo (L. abaculusa, L. abscondita, L. densa, L. dimorpha, L. gibbera, L. hederacea, L. monticola, L. petila, L. rosacea, and L. undulata) (TABLE 2 and FIG. 4). Although there are strong indications that several of the remaining MOTUs could well represent new species as well, at present they are left undescribed, awaiting additional sampling.

Morphology and ecology

A combination of morphological and ecological traits allows a good differentiation of the New Caledonian species. Combinations of morphological, anatomical and ecological characters are graphically represented in FIG. 7.

The Lobophora complex provides an excellent example of the power of molecular-assisted alpha taxonomy (MAAT; Cianciola et al. 2010) in which species are delimited based on molecular data and subsequently the diagnostic value of morphological and ecological characteristics reassessed (see also [START_REF] Verbruggen | Molecular and morphometric data pinpoint species boundaries in Halimeda section rhipsalis (Bryopsidales, Chlorophyta) 1[END_REF][START_REF] Leliaert | DNA-based species delimitation in algae[END_REF]. Even though the current sampling most likely fails dramatically in representing the global species diversity in the genus, several trends with regard to the evolutionary signal of morphological characters stand out. Lineage A composed of five MOTUs, including the newly described species L. rosacea is characterized by a decumbent or fasciculate thallus. Species of lineage B, composed of two species L. nigrescens s.l. and L. australis, is characterized by erect thalli and conspicuous basal holdfasts. Members of the lineage C, including the newly described species L. hederacea, L. undulata, L. monticola and L. abaculusa, are commonly associated with corals and present a predominantly conk-like form.

L. hederacea may also adopt a crustose form especially when found covering specific coral genera (e.g. Seriatopora caliendrum Ehrengerg (1834) and S. hystrix Dana ( 1846)). The species of lineage E, including the species L. dimorpha and L. pachyventera, adopt predominantly a procumbent form. The species of the lineages D and F, including the species L. crassa, L. abscondita, L. gibbera, L. densa, L. asiatica and L. petila, are characterized by a predominant crustose form. L. papenfussii from Bikini Atoll (Marshall Islands) and L. rickeri, from Lord Howe Island (Australia), which presents a crustose form and a thick thallus, may well belong to the lineage F, whose members share the same morphological characteristics (i.e. a crustose form and thick thallus).

Evolutionary perspective and ecological significance of the morphology

The genus Lobophora illustrates the misapprehension of morphological differences for phenotypic plasticity instead of genetic diversity well. Several authors (e.g. De Ruyter van Steveninck et al. 1988;[START_REF] Littler | Caribbean reef plants. An identification guide to the reef plants of the Caribbean, Bahamas, Florida and Gulf of Mexico[END_REF] already observed different growth forms and certainly sensed the existence of different species in relation to the different forms, but nobody ventured to look into this diversity until recently [START_REF] Sun | Taxonomic revision of the genus Lobophora (Dictyotales, Phaeophyceae) based on morphological evidence and analyses rbc L and cox3 gene sequences[END_REF]. The morphological diversity observed within the genus Lobophora was until now considered as the phenotypic plasticity (e.g. [START_REF] Coen | Morphological variation and differential susceptibility to herbivory in the tropical brown alga Lobophora variegata[END_REF]Tanner 1989, De Ruyter van Steveninck et al. 1988;[START_REF] Littler | Caribbean reef plants. An identification guide to the reef plants of the Caribbean, Bahamas, Florida and Gulf of Mexico[END_REF] displayed by a single species, namely Lobophora variegata. Today, three arguments strongly stand against this misconception. First, recent studies including the present one unraveled the hotchpotch of species hidden behind the catch-all species Lobophora variegata. Second, comparison of phylogeny and morphological results revealed the existence of predominant growth forms in each major lineage. Lastly, in a same habitat we may find different species with different forms.

However, one cannot discard phenotypic plasticity off the picture, as we can observe a certain degree of plasticity in every species, with a spectrum of shapes ranging from crustose to erect, but yet again with a predominant form per species. By comparing the morphologies shared by species of a same lineage, we were able to distinguish predominant forms in each lineage. The most basal lineages (A and B) possess predominantly an erect form, the most recent lineages (D-F) present a procumbent to a crustose form, and the intermediate lineage (C) presents a decumbent form. Most likely the ancestral form was a Zonaria-like erect species with a single holdfast, which was also suggested in [START_REF] Sun | Taxonomic revision of the genus Lobophora (Dictyotales, Phaeophyceae) based on morphological evidence and analyses rbc L and cox3 gene sequences[END_REF]. Furthermore, those forms seem to be associated with ecological features. In fact, Lobophora species are found to have a wide variety of habitat and substratum preferences in New Caledonia (e.g. bedrocks, coral rubbles, dead corals and live corals). More remarkably we noticed that this variety of substrata reflected a niche partitioning between the major lineages. For instance, species of lineage B are mostly found growing on sand bottoms, species of lineage C are strongly found in interactions with live corals. These species, present the capacity to bleach and overgrow corals, certainly by the means of secondary metabolites. Species of lineage A are also found in interactions with corals. Species of lineages D to F are mostly found on bedrocks, dead corals or coral rubbles.

Conclusion

The high levels of Lobophora diversity unveiled from a single locality in the Pacific Ocean raises important question with respect to the global diversity of the genus, the distributions and range sizes of the individual species, as well as the mechanisms facilitating co-existence. Current sampling of Lobophora species does not allow to draw far ranging conclusions, but it would appear that individual Lobophora species are restricted to one ocean basin and in this aspect it reminisces the biogeography of the genus Padina, for which there is no or very scanty evidence for species spanning more than one ocean basin. Our analyses included two specimens from the Caribbean Sea, the type locality of L. variegata. Even though the presence of genuine L. variegata in the Indo-Pacific Ocean seems quite unlikely, additional sampling of the Caribbean region is highly needed to precisely determine the identity of L. variegata and assess the species diversity in the Atlantic Ocean. In addition, at present more than half of the MOTUs are recorded only from New Caledonia, but it remains unclear which percentage of the unveiled diversity is really restricted to the study area. An extensive sampling in the Indo-Pacific region is needed to improve our understanding of Lobophora distribution patterns significantly.

C. Vieira is a PhD fellow of the University of Pierre and Marie Curie and Ghent University and is part of MARES, a Joint Doctorate programme selected under Erasmus Mundus coordinated by Ghent University (FPA 2011-0016). Table 1. Comparison of species delimitation analyses. ; f. L. nigrescens s.l..

Tables

Lobophora

  specimens were collected from 41 locations in New Caledonia (FIG. 1). Most of New Caledonia was sampled, except for the remote Entrecasteaux reefs. Sampling sites included the southwest lagoon of Grande Terre (collections between 2004 and 2013), Isle of Pines (BIODIP, November 2005), the Loyalty Islands (BSM-Loyauté, March-April 2005), La Côte Oubliée (CORALCAL1, March 2007), the Chesterfield-Bellona-Bampton area (CORALCAL2, July 2008), Le Grand Lagon Nord (CORALCAL3, February 2009), and different sites along the north west and north east coasts of Grande Terre (CORALCAL4, November-December 2012).

FiguresFigure 1 .

 1 Figures

Figure 2 .

 2 Figure 2. Schematic representation of the various growth forms discerned in Lobophora, with the

Figure 3 .

 3 Figure 3. Schematic representation of a longitudinal and a transverse section of Lobophora,

Figure 4 .

 4 Figure 4. Results of the three species delimitation methods based on the cox3 dataset. Species

Figure 5 .

 5 Figure 5. Boxplots of anatomical variables of New Caledonian Lobophora species; rectangles

Figure 6 .

 6 Figure 6. Lobophora species tree with indication of morphological and ecological characteristics

Figure 7 .

 7 Figure 7. Schematic representation of the ecological (substrate preferences), morphological

Figure

  Figure 8a-l. External morphology of New Caledonian Lobophora species. For new species the

Figure

  Figure 9a-f. Longitudinal (on the left) and transverse (on the right) sections of New Caledonian

Table 2 .

 2 Descriptions of new Lobophora species from New Caledonian.

Acknowledgments

Some samples were collected during several sea campaigns (BIODIP, BSM-Loyauté, CORALCAL1, 2, 3 and 4, MADANG) on board of the IRD research vessel ALIS. We acknowledge the collectors who have contributed to this study: J.-L. Menou, H. Verbruggen, L. Mattio, F. Houlbreque, M. Conord and C. Peignon. We also thank the boat pilots of IRD, Miguel Clarque, Samuel Tereua and Philippe Naudin, for their contribution to this work through their kind and careful assistance.

Supplementary figures

Supplementary tables

Table S1. Origin of the specimens used in this study and their GenBank accession numbers. Table S2. List of primers used in this study.

Table S3. Comparison of species delimitation results for the three methods (GMYC, BI GMYC and ABGD) between the three genes (cox3, rbcL, psbA) on the specimen (left table) and species levels (right table ). Asterisks indicate that the specimens separated by missing sequences are part of the same delimited species. Table S4. Results of the ANOVA of nine anatomical traits for the New Caledonian Lobophora species.

Table S5. Results of the Tukey HSD post-hoc test.

Table S6. Comparison of morphological characters among species of Lobophora.


Pierre-Henri Tournier
email: tournier@ann.jussieu.fr

Frédéric Hecht
email: frederic.hecht@upmc.fr

Myriam Comte
email: comte@ann.jussieu.fr

Finite element model of soil water and nutrient transport with root uptake: explicit geometry and unstructured adaptive meshing

Keywords: finite element method, unstructured mesh adaptation, domain decomposition, plant root uptake

In this paper, we consider a model of soil water and nutrient transport with plant root uptake. The geometry of the plant root system is explicitly taken into account in the soil model. We first describe our modeling approach. Then, we introduce an adaptive mesh refinement procedure enabling us to accurately capture the geometry of the root system and small-scale phenomena in the rhizosphere. Finally, we present a domain decomposition technique for solving the problems arising from the soil model as well as some numerical results.

The paper is organized as follows: section 2 describes our water model. We consider that the root system can be represented as a tree-like network composed of cylindrical root segments. We then define radial and axial water flows on this network that can be coupled to the soil model via a sink term in the Richards equation. The sink term is constructed upon a characteristic function representative of the geometry of the root system. In section 3, the nutrient model is presented in a similar way. Section 4 briefly describes the standard finite element method used in this work. In section 5, the adaptive mesh refinement algorithm is presented: the characteristic function of the root system is computed and then used to construct a metric field in order to drive the mesh adaptation procedure. Since such an approach is computationally intensive, a parallelization technique based on a scalable Schwarz domain decomposition method is used to solve the problems arising from the soil and nutrient models. The procedure is highlighted in section 6.

Throughout this paper the soil domain is denoted by Ω ⊂ R d (d = 2, 3). We consider the evolution of the water potential and nutrient concentration for t ∈ [0, T ], T > 0.

2 The water model

The Richards equation

In soils, water movement is governed by the Richards equation. Richards equation is derived from the continuity equation ∂ θ ∂t

+ ∇.q = S, (2.1)
with θ the volumetric water content, q the macroscopic Darcy flow and S representing sources/sinks. Darcy law relates the water flow to the pressure of the water at any time t:

q = -K∇H, (2.2)
where K is the hydraulic conductivity and H is the total hydraulic head (water potential on weight basis), which can be expressed as

H = h + z. (2.3)
Here, h is the pressure head and comes from a hydrostatic pressure if h > 0 and from a capillary pressure if h < 0. z is the height against the gravitational direction.

The volumetric water content θ and the hydraulic conductivity K are linked to the pressure head h by relationships that depend on the soil properties.

By combining (2.1) and (2.2) we obtain the Richards equation:

∂ t (θ (h)) -∇.(K(h)∇(h + z)) = S in [0, T ] × Ω . (2.4)
Equation (2.4) is subject to the initial condition

h(x, 0) = h 0 (x) in Ω , (2.5)
and the no-flux boundary condition

K(h)∇(h + z).n = 0 on [0, T ] × ∂ Ω , (2.6)
where n denotes the unit outward normal to the boundary of the domain Ω .

The θ (h) and K(h) relationships are given by empirical models whose parameters depend on the soil physical properties. We use the Brooks-Corey model:

Θ (h) := θ (h) -θ m θ M -θ m = h h b -λ := h h b -λ for h ≤ h b 1 for h ≥ h b , K(h) = K s h h b -λ e(λ )
with e(λ

) := 3 + 2 λ , (2.7)
where Θ is the normalized water content.

The parameters of the model are defined as follows:

θ M is the saturated water content.

θ m is the residual water content.

-K s is the saturated hydraulic conductivity.

h b is the bubbling pressure head.

λ is the pore size distribution index.

Experimental evidence has shown that a cycle of wetting-drying of a soil exhibits hysteresis: the water content has different profiles with respect to the wetting and draining processes. This effect can be of importance when considering irrigation or rainfall together with root water uptake. Although hysteresis effects are neglected in the model described here, hysteresis in the soil water retention function θ (h) can be taken into account by including empirical hysteresis models such as [START_REF] Kool | Development and evaluation of closed-form expressions for hysteretic soil hydraulic properties[END_REF] based on main wetting and drying curves.

We introduce the Kirchhoff transformation κ which enables us to reduce the nonlinearity of Richards equation:

κ : h → p h 0 K(p)d p. (2.8)
The new variable p is called the generalized pressure. Previous applications of the Kirchhoff transformation to Richards equation can be found for example in [START_REF] Pop | Error estimates for a time discretization method for the richards' equation[END_REF][START_REF] Varado | Assessment of an efficient numerical solution of the 1d richards' equation on bare soil[END_REF][START_REF] Ji | A generalized transformation approach for simulating steady-state variably-saturated subsurface flow[END_REF][START_REF] Stevens | A scalable and implicit meshless RBF method for the 3D unsteady nonlinear richards equation with single and multi-zone domains[END_REF].

The water content as a function of p is denoted by

M(p) := θ (κ -1 (p)). (2.9)
Using the chain rule, we have

∇p = K(h)∇(h).
(2.10)

Thus the Richards equation reads:

∂ t (M(p)) -∇.(∇p + K(κ -1 (p))∇z) -S = 0 in [0, T ] × Ω . (2.11)
The transformed equation is a semilinear equation in which the nonlinearity in front of the spatial derivative has been eliminated. Using the backward Euler scheme for the time discretization, we are able to write the following weak formulation of the semi-discrete problem: find

p n+1 ∈ H 1 (Ω ) such that ∀v ∈ H 1 (Ω ), Ω M(p n+1 ) -M(p n ) ∆t v + Ω ∇p n+1 ∇v + Ω K(κ -1 (p n+1 ))∇z∇v - Ω Sv = 0. (2.12)
Following the approach suggested in [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF][START_REF] Stevens | A scalable and implicit meshless RBF method for the 3D unsteady nonlinear richards equation with single and multi-zone domains[END_REF], the solution at each time step is obtained iteratively; applying Newton's method to linearize M(p n+1 ) gives the following Newton-like iteration where i is the inner iteration counter for time n + 1:

Ω M ′ (p i )(p i+1 -p i ) + M(p i ) -M(p n ) ∆t v + Ω ∇p i+1 ∇v + Ω K(κ -1 (p i ))∇z∇v - Ω Sv = 0. (2.13)
In practice, as the soil dries the capillary effects get stronger as well as the nonlinearities, while the gravity term becomes of less importance. This allows us to use a simple picard method for the gravity term with no effect on the convergence rate.

The use of the Brooks-Corey model allows us to express the Kirchhoff transformation and its inverse and the transformed functions involved in (2.13) explicitly in a closed form as in [START_REF] Berninger | Domain decomposition methods for elliptic problems with jumping nonlinearities and application to the richards equation[END_REF]:

p = κ(h) = h b -λ e(λ )+1 h h b -λ e(λ )+1 + -λ e(λ )h b -λ e(λ )+1 for h ≤ h b h for h ≥ h b , h = κ -1 (p) =    h b p(-λ e(λ )+1) h b + λ e(λ ) 1 -λ e(λ )+1 for p c < p ≤ h b p for p ≥ h b , M(p) =    θ m + (θ M -θ m ) p(-λ e(λ )+1) h b + λ e(λ ) λ λ e(λ )-1 for p c < p ≤ h b θ M for p ≥ h b , M ′ (p) =    (θ M -θ m ) -λ h b p(-λ e(λ )+1) h b + λ e(λ ) λ λ e(λ )-1 -1 for p c < p ≤ h b 0 for p ≥ h b , K(κ -1 (p)) =    p(-λ e(λ )+1) h b + λ e(λ ) λ e(λ ) λ e(λ )-1 for p c < p ≤ h b 1 for p ≥ h b , (2.14)
Since the transformations (2.14) are obtained analytically, the additional computational cost of employing the Kirchhoff transformation is negligible compared to that of solving the linear system arising from the discretization of (2.13).

Note that the discontinuity of M ′ (p) for p = h b comes from the non-differentiability of the Brooks-Corey function θ (h) at h = h b . However, we do not consider the saturated case h ≥ h b in the examples presented in this paper. Besides, from a numerical point of view, numerical tests have shown that for realistic parameter values the discontinuity is small and does not hinder the convergence of the iterative method.

Root water uptake

Here we consider that the root system is composed of cylindrical root segments. It can then be represented as a series of interconnected nodes forming a network of segments Σ , each segment with its own parameters (radius, conductivity, etc.). Such a representation can be generated by the Matlab code RootBox [START_REF] Leitner | A dynamic root system growth model based on l-systems[END_REF], which implements a root system growth model based on L-Systems. RootBox is a root architectural model that explicitly simulates the architecture of root systems in the 3D space, using a set of growth rules which are applied to a series of root types or classes, with each root type having its own characteristic set of growth parameters such as root elongation rate or branching density. The algorithm computes elongation and branching of the roots according to the initial growth speed, lengths of apical and basal zones as well as internodal distances, maximal number of branches and branching angles. Growth direction can follow different types of user defined tropisms. The model has a stochastic component in that all parameters can be given with mean and standard deviation. Fig. 1 shows an example of a maize root system generated by RootBox.

In the following we define water flows on the root network and describe the coupling with the soil model as was done in [START_REF] Doussan | Water uptake by plant roots: II -modelling of water transfer in the soil root-system with explicit account of flow within the root system -comparison with experiments[END_REF].

We make the same assumptions as in [START_REF] Doussan | Water uptake by plant roots: II -modelling of water transfer in the soil root-system with explicit account of flow within the root system -comparison with experiments[END_REF], which are explained in [START_REF] Doussan | Modelling of the hydraulic architecture of root systems: An integrated approach to water absorption -model description[END_REF]. The model describing root water uptake can be found in [START_REF] Fiscus | The interaction between osmotic-and pressure-induced water flow in plant roots[END_REF][START_REF] Landsberg | Water movement through plant roots[END_REF]: water flow in roots and between soil and roots is described by the transpiration-cohesiontension mechanism and follows an Ohm's law analogy. Let us recall the hypotheses made in [START_REF] Doussan | Modelling of the hydraulic architecture of root systems: An integrated approach to water absorption -model description[END_REF]: First, the influence of solutes on flow is neglected, because during periods of active transpiration, the hydrostatic pressure gradient rather than the osmotic potential gradient is the effective driving force for flow. The second hypothesis consists in neglecting the capacitive effect of the roots and considering only steady-state flow, because water stored in roots is generally small compared to transpiration requirements. Thus, for a cylindrical root segment of radius r and length l and following [START_REF] Landsberg | Water movement through plant roots[END_REF], we can define the volumetric radial water flow into the root from the soil J r and the longitudinal flow up the root in the xylem J x as Although these simplifications are made in our model as well, the model can be extended by taking into account osmotic gradients and capacitive effects of roots.

J r = L r s r (h s -h r ), J x = -K r dh r dl , (2.15)
Equations (2.15) giving radial and longitudinal flows can be used to formulate a water mass balance equation for a given root node i of parent node p in the tree-like structure as depicted in Fig. 2:

J x,i = ∑ j∈childs(i) J x, j + J r,i , (2.16)
which can be written as

-K r,i h r,p -h r,i l i = -∑ j∈childs(i) K r, j h r,i -h r, j l j + L r,i 2πr i l i (h s,i -h r,i ) + (h s,p -h r,p ) 2 .
(2.17)

Here K r,i , L r,i , r i and l i refer to the root segment (p, i) while K r, j and l j relate to the root segment (i, j). h s,i and h r,i are the soil water potential at root node i and the xylem water potential at root node i respectively. We approximate the potentials h s and h r for segment (i, p) by averaging their value at the two nodes i and p. Parameters L r and K r are given for each segment and can depend on various data such as root type and age. Writing (2.17) for every node in the tree-like structure, the xylem water potential vector (h r,i ) i is then solution of a linear system, with the right-hand side containing the soil factors represented by the h s,i . At the root collar, we can prescribe the transpiration flow or the xylem potential with a Neumann or Dirichlet boundary condition respectively. We can follow the same approach as in [START_REF] Doussan | Water uptake by plant roots: II -modelling of water transfer in the soil root-system with explicit account of flow within the root system -comparison with experiments[END_REF][START_REF] Javaux | Use of a three-dimensional detailed modeling approach for predicting root water uptake[END_REF]: In the case of a flux-type boundary condition, stress may occur when the evaporative demand cannot be met by the soil. In such a case, a maximum Other models could also be considered, such as [START_REF] Tuzet | A coupled model of stomatal conductance, photosynthesis and transpiration[END_REF] where stomatal response to a drying soil is modeled by a logistic function with empirically determined parameters.

p i j1 j2 J x, j 1 J x, j 2 J x,i J r,i
In order to take the radial water uptake flows into account in the soil water model, a sink term S in the Richards equation is defined in the domain. Since the sink term represents root uptake flows in the 3D (or 2D) space, we construct S through a characteristic function of the root system f c representative of its geometry using the distance function to the root network Σ . The characteristic function f c can be seen as a smooth approximation of the 1D root network Σ , taking the values 1 at the root and 0 away from the root, with a smooth change of width ε in between. The purpose of the characteristic function is threefold: define a regularization of the delta function representing the network of segments Σ , construct a sink term matching the volume occupied by the roots by using the diameter of the root as the width of the regularization, and drive the adaptive mesh refinement procedure.

The function f c representative of the geometry of the root system in the domain is constructed as follows:

-For a point x of the domain Ω the distance d from x to the root is computed:

d(x) = min s∈Σ d s (x), (2.18)
with Σ the set of root segments in the tree-like network. For each root segment s, the distance d s (x) from the point x to the segment s is easily computed using distance from line and point routines.

-The distance function d is then used to compute the characteristic function. There is a variety of admissible transformations that we can use, and we choose the following:

f c (x) = f d (d(x)) = 1 -tanh 6d(x) ε . (2.19)
We can choose ε to be equal to the radius of the root.

We can now build the sink term in the Richards equation. In order to ensure that the sink term in the soil model corresponds to the volumetric radial flow in the network model, we need to introduce a scale factor depending on the choice of f d .

Let us consider the case of a cylindrical root segment s, formed by the nodes i and j. The corresponding radial flow is

J r = L r 2πr r l r (h s,i -h r,i ) + (h s, j -h r, j ) 2 .
(2.20)

We want the integral of the corresponding sink term S over the domain to be equal to the outflow rate, i.e.

Ω S = -J r .

(2.21)

If f c is the characteristic function of the single root segment as defined above, using cylindrical coordinates we get (in the 3D case)

Ω f c ≃ 2πl r R 0 r f d (r)dr (2.22)
with R >> ε. The approximation error coming from the truncature in the integral is negligible for usual choices of f d . Let us define Tr as

T r = R 0 r f d (r)dr. (2.23)
We then define the sink term S as

S = -f c L r r r T r h l , (2.24)
where h l only depends on the longitudinal coordinate and linearly interpolates h sh r along the segment. Thus, we have

Ω S = - Ω f c L r r r T r h l = -L r 2πr r l r (h s,i -h r,i ) + (h s, j -h r, j ) 2 = -J r .
(2.25)

The extension to the whole root system is straightforward.

Since the characteristic function does not correspond to an arrangement of perfect cylindrical root segments due to its shape at root tips or in-between root segments, a modified approach consists in adjusting the surface areas of the root segments in the definition of the radial uptake flows in the root network model so that for each root segment, the volumetric uptake flow is equal to the actual contribution of the segment to the global sink term in the soil model. This approach ensures that the amount of water depleted in the soil water model is equal to the transpiration rate in the network model, although the difference is minimal in actual computations.

The coupling between the tree-like model and the soil water model consists in iteratively solving the two problems until convergence. Let h t i s be the soil matric potential distribution at time t i , h k s and h k r the soil and xylem matric potentials at inner iteration k and time t i+1 . The coupling algorithm reads as follows:

1. h 0 s = h t i s . 2. Solve the linear system of the tree-like model derived from (2.17) with soil factors h k s , obtain h k r on the root network.

3. Compute the sink term S as in (2.24) using h k s and h k r . 4. Perform an inner iteration of (2.13), obtain h s in the soil domain.

h k+1

s = h k s + α k (h s -h k s )
, where α k is an under-relaxation parameter that ensures convergence of the system. 6. If ||h sh k s || > τ, go to 2. with k ← k + 1. Fig. 3 gives an overview of the water model through an example.

The nutrient model

Here we consider the evolution of the concentration c of a nutrient ion N in the soil solution, governed by diffusion, mass flow, adsorption in the soil solid phase and root uptake.

The convection-diffusion equation

The convection-diffusion equation expresses the nutrient mass balance and can be written in its conservative form: -A is the diffusion coefficient of N in the soil. For simplicity we consider a modified diffusion coefficient with a tortuosity factor in place of the effective dispersion coefficient tensor that is usually employed for flow in porous media, although the model can easily be extended to account for dispersibility since water flow is explicitly considered. Following [START_REF] Barraclough | The determination of ionic diffusion coefficients in field soils. I. diffusion coefficients in sieved soils in relation to water content and bulk density[END_REF], A is given by

∂ t (θ c + ϕ(c)) = ∇. A∇c -qc) + S c (c) in [0, T ] × Ω , ( 3
A = A 0 θ f l , (3.2)
with A 0 the diffusion coefficient of N in free water and f l the so-called tortuosity factor:

f l = f 1 θ + f 2 for θ ≥ θ l , θ ( f 1 θ l + f 2 ) θ l for θ < θ l , (3.3)
where f 1 , f 2 and θ l are parameters depending on soil properties. -ϕ is an adsorption/desorption isotherm relating the amount of adsorbed N to the equilibrium concentration of N in solution; we use the Freundlich adsorption isotherm [START_REF] Mcgechan | Sorption of phosphorus by soil, part 1: Principles, equations and models[END_REF], defined by Here θ + ϕ ′ (c) is the buffer power and represents the ability of the soil to resupply nutrients as plants take it up from the soil solution.

ϕ(c) = κc b , (3.4
The Lagrangian form is derived by dividing by θ + ϕ ′ (c) and defining the material derivative

Dc Dt = ∂ t c + q θ + ϕ ′ (c) .∇c. (3.8)
We can use the method of characteristics to handle the convective part. The velocity field is a = q θ +ϕ ′ (c) . Following [START_REF] Saaltink | Mass balance errors when solving the convective form of the transport equation in transient flow problems[END_REF], we can use the approach yielding (3.7) at the semi-discrete level to reduce mass balance errors. Using the implicit Euler scheme for the time discretization of (3.1) gives

θ n+1 c n+1 + ϕ(c n+1 ) -θ n c n -ϕ(c n ) ∆t = ∇. A n+1 ∇c n+1 -q n+1 c n+1 + S c (c n+1 ). (3.9)
Our implicit time discretization of Richards equation (2.12) corresponds to

θ n+1 -θ n ∆t + ∇.q n+1 = S n+1 . (3.10)
Multiplying (3.10) by c n+1 and substracting it from (3.9), we are able to write the weak formulation of the semidiscrete problem: find

c n+1 ∈ H 1 (Ω ) such that ∀v ∈ H 1 (Ω ), Ω θ n c n+1 + ϕ(c n+1 ) -(θ n c n ) • X -ϕ(c n • X) ∆t v + Ω A n+1 ∇c n+1 ∇v + Ω S n+1 c n+1 v - Ω S c (c n+1 )v = 0 (3.11)
with the approximated characteristics X = x -∆ta. Applying Newton's method and denoting by i the inner iteration counter for time n + 1, we get

Ω θ n c i+1 + ϕ(c i ) + ϕ ′ (c i )(c i+1 -c i ) -(θ n c n ) • X -ϕ(c n • X) ∆t v + Ω A n+1 ∇c i+1 ∇v + Ω S n+1 c i+1 v - Ω S c (c i )v - Ω S ′ c (c i )(c i+1 -c i )v = 0.
(3.12)

Nutrient uptake

The root acts as a selective membrane for ion uptake. Nutrient uptake by roots is given by a model of enzyme kinetics, relating in this case the root uptake rate of N to its concentration at the root surface. Here we use the following Michaelis-Menten model:

h(c) = F m c K m + c , (3.13)
where h(c) is the uptake rate and F m , K m > 0 are parameters of the model. This model of active nutrient uptake is taken from [START_REF] Barber | Soil nutrient bioavailability: a mechanistic approach[END_REF]. The hypothesis that active uptake becomes more important under low nutrient supply, while the transpiration driven mass flow dominates for higher concentrations, is proposed in [START_REF] Šimůnek | Modeling compensated root water and nutrient uptake[END_REF]. For a cylindrical root segment of radius r and length l and assuming that (3.13) we can define the volumetric nutrient uptake rate by

J N = 2πrlh(c s ), (3.14)
where c s is the concentration of N at the root surface.

Similarly to the sink term (2.24) in the Richards equation, the sink term S c representing nutrient uptake in the convection-diffusion equation is constructed as follows:

S c = -f c r T r h(c). (3.15)
The model can easily be adapted to implement other nutrient uptake models. For example, in [START_REF] Somma | Transient three-dimensional modeling of soil water and solute transport with simultaneous root growth, root water and nutrient uptake[END_REF][START_REF] Hopmans | Current capabilities and future needs of root water and nutrient uptake modeling[END_REF] the solute uptake term is defined as

S ′ (c) = εSc + (1 -ε)S c (c),
where ε ∈ [0, 1] is a coefficient partitioning total uptake between passive uptake Sc where solute enters the root dissolved in water and active uptake S c (c), which could also be described following for example [START_REF] Kochian | Potassium transport in corn roots: I. resolution of kinetics into a saturable and linear component[END_REF] by a Michaelis-Menten-type kinetic with a linear component.

Finite element formulation

In this section we describe briefly the Galerkin P1 finite element approximation of problems (2.13) and (3.12). Let T h be a mesh of the domain Ω . Let

V h = u h ∈ H 1 (Ω ) u h | K ∈ P 1 , ∀K ∈ T h . (4.1)
Spatial discretization of Equation (2.13) leads to the following discrete variational problem: find

p h ∈ V h such that ∀v h ∈ V h , M ′ (p i h )p h , v h + ∆t (∇p h , ∇v h ) = M ′ (p i h )p i h , v h -M(p i h ), v h + (M(p n h ), v h ) -K(κ -1 (p i h ))∇z, ∇v h + S i h , v h . (4.2)
Spatial discretization of Equation (3.12) leads to the following discrete variational problem: find

c h ∈ V h such that ∀v h ∈ V h , (θ n h c h , v h ) + ϕ ′ (c i h )c h , v h + A n+1 h ∇c h , ∇v h + ∆t S n+1 h c h , v h -∆t S ′ c (c i ) h c h , v h = -ϕ(c i h ), v h + ϕ ′ (c i h )c i h , v h + ((θ n h c n h ) • X, v h ) + (ϕ(c n h • X), v h ) + ∆t S c (c i ) h , v h -∆t S ′ c (c i ) h c i h , v h . (4.3)
Numerical resolution of systems (4.2) and (4.3) is carried out using the finite element software FreeFem++ [START_REF] Hecht | New development in freefem++[END_REF].

Mesh adaptation

The soil domain Ω is first represented by a regular initial simplicial finite element mesh. Since the characteristic function of the root system f c is poorly represented on the initial mesh, we refine it iteratively using anisotropic mesh adaptation. Since we expect high gradients and small-scale phenomena to be localized near the roots (i.e. where f c exhibits strong variations), this type of a priori refinement is adequate. The main steps of the adaptive procedure are as follows: First, we compute f c for each node of the mesh. Then we define a nodal based anisotropic metric from the Hessian of the function f c . Finally, the mesh is adapted using the size and streching of elements provided by the metric. This procedure is repeated iteratively. In 2D, we use the built-in adaptive remesher of FreeFem++. For 3D simulations, FreeFem++ is interfaced with mshmet for computing the Hessian-based anisotropic metric and with the anisotropic fully tetrahedral automatic remesher Mmg3d [START_REF] Dobrzynski | MMG3D: User Guide[END_REF] which uses anisotropic Delaunay kernel and local mesh modifications based on a combination of edge flips, edge collapsing, node relocation and vertex insertion operations to adapt the mesh. Fig. 4 illustrates the mesh adaptation process in a 2D simulation.

Domain decomposition

Sinces meshes generated as described in section 5 require a considerable number of nodes to be able to adequately resolve the geometry of complex root systems, linear systems resulting from the discrete problems (4.2) and (4.3) can be quite large. In order to reduce computation time, we opted for a parallel divide-and-conquer technique with an additive Schwartz overlapping domain decomposition method. The initial computational domain is partitioned by metis [START_REF] Karypis | A fast and high quality multilevel scheme for partitioning irregular graphs[END_REF] into a number of subdomains, on which local variational problems are defined. A two-level coarse grid preconditioner taken from [START_REF] Jolivet | High performance domain decomposition methods on massively parallel architectures with freefem++[END_REF] is also used to improve the convergence of the domain decomposition method. Numerical tests are conducted in 2D and in 3D in order to assess the efficiency of the two-level preconditioner compared to a classical one-level preconditioner. We consider one inner iteration of (4.2). In the 2D case, the mesh is composed of 204 331 vertices and 407 839 triangles. The 3D mesh is composed of There are several iterative algorithms used to obtain the numerical solution. In the water model, the outer loop consists in a fixed point algorithm solving alternatively the root problem (2.18) and the soil problem (4.2). Solving the linear system resulting from the discrete soil problem (4.2) using the domain decomposition method presented in this section constitutes an inner loop. Local problems defined on each subdomain are also solved iteratively using a conjugate gradient method, and thus the complete algorithm consists in three nested loops. We can take advantage of the iterative nature of the domain decomposition and linear solvers by using adaptive stopping criteria in order to further reduce the computational time. Here we use simple heuristics expressing that there is no need to continue with iterations in the inner loop once the error from the outer loop starts to dominate. More elaborate stopping criteria can be used, see for example [START_REF] Vohralík | A posteriori error estimates, stopping criteria, and adaptivity for two-phase flows[END_REF] where adaptive stopping criteria based on a posteriori error estimates are derived. the initial water potential in the soil domain is in hydrostatic equilibrium: h 0 = -15 m -z.

For simplicity, all root parameters are taken constant across the whole root system. Numerical values of L r and K r for maize are taken from [START_REF] Doussan | Modelling of the hydraulic architecture of root systems: An integrated approach to water absorption -model description[END_REF]:

-Root radius is set to 5.0 × 10 -4 m.

-L r = 1.923 08 × 10 -4 d -1 , K r = 4.32 × 10 -8 m 3 d -1 .

As an illustration of the nutrient model, we consider the transport and uptake of nitrate. Parameters for equation (3.2) are taken from [START_REF] Barraclough | The determination of ionic diffusion coefficients in field soils. I. diffusion coefficients in sieved soils in relation to water content and bulk density[END_REF], Michaelis-Menten constants for maize are taken from [START_REF] Barber | Soil nutrient bioavailability: a mechanistic approach[END_REF]:

-A 0 = 1.6416 × 10 -4 m 2 d -1 , f 1 = 1.58 , f 2 = -0.17, θ l = 0.12. -F m = 8.64 × 10 -3 mol m -2 d -1 , K m = 2.5 × 10 -2 mol m -3 .
the homogeneous initial concentration of nitrate in the soil solution is c 0 = 5 mol m -3 .

We consider that adsorption of nitrate in the soil solid phase is negligible: ϕ = 0.

The first numerical simulation involves the 20-days-old maize root system generated by the Matlab code Root-Box depicted in Fig. 1. The soil domain is of dimensions 0.4 m × 0.4 m × 0.4 m. No-flux boundary conditions are imposed on the boundaries of the soil domain. A constant transpiration rate equal to 1.44 × 10 -4 m 3 d -1 is imposed at the root collar. The time step ∆t is taken constant equal to 0.05 d. Numerical results are depicted in Figs. 5 and6. Fig. 5 shows the Darcy flux q in the vicinity of the roots. In the beginning of the simulation, root water uptake is still relatively evenly distributed over the dense upper portion of the root system (left picture, t = 4.8h). As time passes, the uptake pattern is modified. The soil dries in the dense root zone and the root system takes up water from wetter zones (mostly in the deeper part of the soil profile) in order to maintain a constant transpiration rate (right picture, t = 7.3d). Fig. 6 shows high gradients developing in the vicinity of the roots as the soil dries.

The second example illustrates how we are able to integrate root growth and chemotropism in 2D by coupling the model with the implementation of growth and tropisms in RootBox. RootBox can simulate root tip response to mechanical soil heterogeneities as well as various types of tropisms like gravitropism, hydrotropism or chemotropism. The specific growth behaviour can be chosen for every root type. In RootBox, the implementation of tropisms consists in computing the new growth direction by random minimization of an objective function: for each active root tip, several rotations are randomly computed and the one that leads to minimizing the objective function is chosen. Different types of tropisms are realized by choosing appropriate objective functions depending on soil properties (water content, nutrient concentration). In this example, a combination of gravitropism and chemotropism is achieved by defining the objective function as -λ c + z where λ > 0 represents the relative strength of chemotropism. At each time step, an iteration of the following coupling algorithm is performed:

1. RootBox finds the best growth direction for each active root tip through multiple evaluations of the objective function depending on the spatial concentration c. RootBox is interfaced with the FreeFem++ finite element code so that the values of the concentration can be determined by interpolation on the mesh. Then, a new tree-like network with new root segments is obtained. 2. A mesh adapted to the new characteristic function f c is obtained from the previous mesh by the mesh adaptation procedure described in section 5. 3. Finite element functions (namely the current soil water potential and N concentration) are interpolated from the previous mesh to the new mesh. 4. Solve (2.12) and (3.11) and obtain the new soil water potential and N concentration distributions. Fig. 7 depicts some results of such a simulation with the initial nitrate concentration set to a linear profile varying from 0 mol m -3 at the top left corner of the domain to 10 mol m -3 at the bottom right corner. Notice the accumulation of nitrate around some of the roots in the bottom right: as the soil dries out, radial soil-root water flow increases in the bottom right where the soil is wetter, resulting in the mass flow of nitrate bringing more than the root can take up.

Conclusion

In this paper we presented a model of soil water and nutrient transport with plant root uptake. A characteristic function of the geometry of the root system was used to construct accurate sink terms corresponding to water and nutrient uptake by roots. An emphasis was put on the spatial discretization with an adaptive mesh refinement procedure producing meshes that are able to resolve the complex geometry of the root system together with smallscale phenomena occuring in the rhizosphere. A parallel finite element method was then presented using a twolevel Schwarz domain decomposition method to solve the potentially large systems arising from the discretization. Numerical experiments were conducted in two and three spatial dimensions to illustrate the capabilities of the model. Future work will consist in using the diffuse domain approach to approximate the actual surface of the roots in the volume mesh without relying on a surface mesh. The level set representing the root surface splits the domain Ω into two subdomains Ω s and Ω r . Following a monolithic approach, different problems are defined on the same mesh in the soil domain Ω s and in the root system Ω r using a phase field variable. This approach allows us to make use of the same adaptive volume remeshing procedure that was described in this paper.

Fig. 1

 1 Fig. 1 Example of a 20-days-old maize root system generated by RootBox composed of 10611 segments

Fig. 2

 2 Fig.2Water mass balance for root node i

Fig. 3

 3 Fig. 3 Overview of the water model

  ) where κ > 0 and b ∈ (0, 1) are fitting parameters of the model, -S c (c) represents sources/sinks. Equation (3.1) is subject to the initial condition c(x, 0) = c 0 (x) in Ω , (3.5) and the no-flux boundary condition A∇cqc).n = 0 on [0, T ] × ∂ Ω . (3.6) The convective form of equation (3.1) is obtained by multiplying equation (2.1) by c and substracting it from equation (3.1). This gives (θ + ϕ ′ (c))∂ t c = ∇. A∇c)q.∇c -Sc + S c (c). (3.7)

Fig. 4

 4 Fig. 4 2D example of the mesh adaptation process: representation of the function f c (left) defined on the adapted mesh (right).

7

 7 Numerical resolutionThe purpose of the numerical examples presented in this section is to illustrate the capabilities of the numerical model.Numerical values used in the examples are as follows:parameters for a clay soil are θ m = 0.068, θ M = 0.38, λ = 0.17, h b = -0.4 m, K s = 0.144 m d -1 .

Fig. 5

 5 Fig. 5 Soil water flow near the root system of a 20-days-old maize plant

Fig. 7

 7 Fig. 7 2D simulation example with root growth and chemotropism: snapshots of nitrate concentration at different time steps. In white, isosurface 0.5 of the characteristic function f c .

Acknowledgements We thank Pierre Jolivet for his valuable insights on FreeFem++ and domain decomposition methods.