Mélodie Ollivier

Vincent Lesieur

Johannes Tavoillot

Fanny Bénetière

Jean-François Martin

Marie-Stéphane Tixier

Metabarcoding and ecological interaction networks for selecting candidate biological control agents

Keywords: High-throughput sequencing, host range, food web, weed biocontrol, Common Sowthistle

come   L'archive ouverte pluridisciplinaire

2 Abstract 1. Classical biological control can be used to decrease the density of invasive species to below an acceptable ecological and economic threshold. Natural enemies specific to the invasive species are selected from its native range and released into the invaded range. This approach has drawbacks, despite the performance of specificity tests to ensure its safety, because the fundamental host range defined under controlled conditions does not represent the actual host range in natura, and these tests omit indirect interactions within community.

2. We focus on Sonchus oleraceus (Asteraceae), a weed species originating from Western Palearctic that is invasive worldwide and notably in Australia. We explore how analyses of interaction network within its native range can be used to 1) inventory herbivores associated to the target plant, 2) characterize their ecological host ranges, and 3) guide the selection of candidate biocontrol agents considering interactions with species from higher trophic levels. Arthropods were collected from plant community sympatric to S. oleraceus, in three bioclimatic regions, and interactions were inferred by a combination of molecular and morphological approaches.

3. The networks reconstructed were structured in several trophic levels from basal species (diversified plant community), to intermediate and top species (herbivorous arthropods and their natural enemies). The subnetwork centered on S. oleraceus related interactions contained 116 taxa and 213 interactions. We identified 47 herbivores feeding on S. oleraceus, 15 of which were specific to the target species (i.e. Generality value equal to 1). Some discrepancies with respect to published findings or conventional specificity tests suggested possible insufficient sampling effort for the recording of interactions or the Introduction Biological invasions are currently threatening biodiversity to an unprecedented extent [START_REF] Bellard | Global patterns in threats to vertebrates by biological invasions[END_REF][START_REF] Seebens | Global trade will accelerate plant invasions in emerging economies under climate change[END_REF][START_REF] Vitousek | Human Domination of Earth's Ecosystems[END_REF]. When invasive species disrupt the ecological or economic balance, action is required to control their negative impact. Chemical control methods are widely used in such situations, but classical biological control (CBC) constitutes a possible alternative. CBC involves the release of natural enemies, specific to the target organism and originating from its native range, to keep the density of the invasive species below an economically and ecologically acceptable threshold [START_REF] Keane | Exotic plant invasions and the enemy release hypothesis[END_REF][START_REF] Mcfadyen | Biological control of weeds[END_REF][START_REF] Van Driesche | Classical biological control for the protection of natural ecosystems[END_REF]. CBC is considered more sustainable than chemical control [START_REF] Peterson | Comparing biocontrol and herbicide for managing an invasive non-native plant species: Efficacy, non-target effects and secondary invasion[END_REF], although the introduction of biocontrol agents (BCA) into a new territory may itself represents a risk for the recipient communities [START_REF] Barratt | The status of biological control and recommendations for improving uptake for the future[END_REF][START_REF] Hinz | How Safe Is Weed Biological Control? A Global Review of Direct Nontarget Attack[END_REF][START_REF] Suckling | What Magnitude Are Observed Non-Target Impacts from Weed Biocontrol?[END_REF]. Once introduced, the BCA may affect non-target species, especially if it lacks specificity [START_REF] Müller-Schärer | Classical biological control: exploiting enemy escape to manage plant invasions[END_REF]. Assessing the host range of a candidate BCA is, thus, crucial, to anticipate such risks. Most of the host-specificity tests performed to assess this risk are conducted under standardized conditions, through choice/nochoice experiments over a range of targets selected according to the centrifugal phylogeny approach [START_REF] Briese | Translating host-specificity test results into the real world: The need to harmonize the yin and yang of current testing procedures[END_REF][START_REF] Wapshere | Recent developments in biological control of weeds[END_REF].

Recent reviews recognized the success of such experiments for limiting the undesirable unintentional effects of the CBC of weeds [START_REF] Hinz | How Safe Is Weed Biological Control? A Global Review of Direct Nontarget Attack[END_REF][START_REF] Hinz | A global review of target impact and direct nontarget effects of classical weed biological control[END_REF]. However, the cumbersome nature of these tests reduces the range of species that can be screened. The candidate BCA are selected through preliminary field monitoring that may miss a species of interest. Furthermore, as the fundamental host range of a species (defined under controlled conditions) is thought to be broader than the host range actually observed in the field (known as the realized host range; [START_REF] Louda | Nontarget effects--the Achilles' heel of biological control? Retrospective analyses to reduce risk associated with biocontrol introductions[END_REF][START_REF] Schoonhoven | Insect-plant biology: from physiology to evolution[END_REF][START_REF] Sheppard | Advances in molecular ecology: tracking trophic links through predator-prey food-webs[END_REF], these tests tend to overestimate the risk and lead to the rejection of candidate BCA based on interactions that would not occur in the field (false positives) (e.g. [START_REF] Groenteman | St. John's wort beetles would not have been introduced to New Zealand now: A retrospective host range test of New Zealand's most successful weed biocontrol agents[END_REF]. Most CBC programs use specificity tests under controlled conditions as proxies for field conditions due to the complexity of trophic interaction assessments in the field but this leaves room for improvement.

The characterization of ecological interactions among communities of plants and arthropods in natura is challenging, as it traditionally requires direct observations, the rearing of specimens and considerable taxonomic expertise, rendering the process impractical for large-scale studies.

Recent advances in molecular approaches, such as the combination of DNA metabarcoding on gut content or feces and high-throughput next-generation sequencing (NGS), have opened up new opportunities to track the host range of arthropods in natura with both a high taxonomic resolution and high sensitivity [START_REF] Derocles | Chapter One -Biomonitoring for the 21st Century: Integrating Next-Generation Sequencing Into Ecological Network Analysis[END_REF][START_REF] Frei | Molecular analysis indicates high levels of carabid weed seed consumption in cereal fields across Central Europe[END_REF][START_REF] Wirta | Complementary molecular information changes our perception of food web structure[END_REF][START_REF] Zhu | Seeing is believing? Comparing plant-herbivore networks constructed by field co-occurrence and DNA barcoding methods for gaining insights into network structures[END_REF]. Even interactions that are very difficult to observe, such as host-parasitoid associations, can be detected by such methods [START_REF] Gariepy | A molecular diagnostic tool for the preliminary assessment of host-parasitoid associations in biological control programmes for a new invasive pest[END_REF][START_REF] Hrček | What do molecular methods bring to host-parasitoid food webs?[END_REF].

This approach can be used to reconstruct networks of trophic interactions directly from studies in the field, and provides an analytical framework particularly relevant to studies of complex species assemblages. Network ecology do not only depicts species interactions, but provides elements for the understanding of recurrent patterns of antagonistic interactions between plants and herbivores, such as specialization or compartmentation [START_REF] Lewinsohn | Structure in plant-animal interaction assemblages[END_REF][START_REF] Thébault | Stability of Ecological Communities and the Architecture of Mutualistic and Trophic Networks[END_REF]. In CBC against invasive weeds, analyses of ecological networks have been used to assess the extent to which a BCA fits into a recipient community. Such methods provided a way to quantify the direct impact of biological control on non-target plants [START_REF] Memmott | Food webs as a tool for studying nontarget effects in biological control[END_REF], and its indirect impact on other species at higher trophic levels [START_REF] Carvalheiro | Apparent competition can compromise the safety of highly specific biocontrol agents[END_REF][START_REF] Louda | Ecological Effects of an Insect Introduced for the Biological Control of Weeds[END_REF][START_REF] Pearson | Indirect effects of host-specific biological control agents[END_REF]. Such studies have highlighted the usefulness of network ecology for evaluating the impact of BCA after their introduction [START_REF] Memmott | Food webs: a ladder for picking strawberries or a practical tool for practical problems?[END_REF][START_REF] Willis | The potential for indirect effects between a weed, one of its biocontrol agents and native herbivores: A food web approach[END_REF], but interaction network analysis can also be used for the upstream assessment of potential candidate BCA, in a more systematic process [START_REF] Ollivier | Characterizing ecological interaction networks to support risk assessment in classical biological control of weeds[END_REF]. Adopting a network approach as prerequisite step of the CBC program, can provide a rapid screening of the ecological host range of potential agents to be tested in priority. This can also inform about species functional properties through the position and connexions the species have in the network, independently of its taxonomic assignation, which would confer a strong predictive power of the interactions possibly occurring in a novel bioclimatic region [START_REF] Todd | Using qualitative food webs to predict species at risk of indirect effects from a proposed biological control agent[END_REF].

Indeed, the choice of BCA should also take into account indirect effects on the recipient community due to interactions with higher trophic levels in the network, i.e. natural enemies [START_REF] Hinz | How Safe Is Weed Biological Control? A Global Review of Direct Nontarget Attack[END_REF][START_REF] Memmott | Food webs as a tool for studying nontarget effects in biological control[END_REF]. If comparable enemies than those identified in the native range are present in the range of introduction, new interactions might be created with BCA, resulting in disturbances in the ecological network through indirect interactions, e.g. apparent competition [START_REF] Carvalheiro | Apparent competition can compromise the safety of highly specific biocontrol agents[END_REF][START_REF] López-Núñez | Fourtrophic level food webs reveal the cascading impacts of an invasive plant targeted for biocontrol[END_REF].

The objective of this study was to determine how the analysis of interaction networks could be used to support the selection of candidate BCA for the common sowthistle, Sonchus oleraceus L.

(Asteraceae). This plant is native to Western Europe and Northern Africa [START_REF] Boulos | Revision systematique du genre Sonchus L. s.l. VI. Sous-genre. 3. origosonchus. genres Embergeria, Babcockia et -Taeckholmia. species exclusae et dubiae[END_REF][START_REF] Hutchinson | THE BIOLOGY OF CANADIAN WEEDS.: 63. Sonchus asper (L.) Hill and S. oleraceus L. Can[END_REF] and is the most widely naturalized terrestrial plant worldwide [START_REF] Pyšek | Naturalized alien flora of the world: species diversity, taxonomic and phylogenetic patterns, geographic distribution and global hotspots of plant invasion[END_REF]. In Australia, it has become a weed of major concern in cropping systems [START_REF] Llewellyn | Impact of weeds on Australian grain production : the cost of weeds to Australian grain growers and the adoption of weed management and tillage practices[END_REF][START_REF] White | Germination, emergence, and persistence of Sonchus oleraceus, a major crop weed in subtropical Australia[END_REF]. Aside the development of resistance to multiple herbicides [START_REF] Adkins | Weeds resistant to chlorsulfuron and atrazine from the north-east grain region of Australia[END_REF][START_REF] Jalaludin | Glyphosate and 2,4-D amine resistance in common sowthistle (Sonchus oleraceus) and fleabane (Conyza bonariensis) in the northern grain growing region of Australia[END_REF][START_REF] Meulen | Survey of glyphosate resistance in common sowthistle (Sonchus oleraceus) across the Australian Northern Grains Region[END_REF], the control of this weed is complex as it is extremely prolific and seeds can germinate all year round when sufficient rainfalls occurs. Sonchus oleraceus rapidly dominates crops, reducing yield and contaminating harvested grain [START_REF] Llewellyn | Impact of weeds on Australian grain production : the cost of weeds to Australian grain growers and the adoption of weed management and tillage practices[END_REF]. A CBC program was therefore initiated in 2017, to identify candidate BCA. In this context, an analysis of ecological networks, based on direct field observations and high-throughput DNA metabarcoding, was performed. Our objectives were to 1) establish an inventory of arthropods feeding on S. oleraceus, and assess the contribution of the approach relatively to classical procedures, 2) delineate the ecological host range for herbivores feeding on S. oleraceus and point out candidate BCA, and 3) identify the trophic interactions of the candidate BCA with natural enemies, and consider their implications for the CBC program.

Materials and Methods

Sampling design

We maximized the species diversity and associated interactions, through a maximum variation design with three bioclimatic regions in France (semi-oceanic, Mediterranean and continental climates) [START_REF] Ceglar | Observed Northward Migration of Agro-Climate Zones in Europe Will Further Accelerate Under Climate Change[END_REF] and three successive sampling dates (April, May and June 2018).

Sampling was carried out from 10 A.M. to 4 P.M., by varying climatic conditions (wet, cloudy to sunny weathers and temperatures ranged between 10°C and 29°C). These variations did not affect our ability to capture arthropods. For each bioclimatic region and date, we employed an opportunistic sampling strategy to collect plants from several ruderal and agricultural sites, covering the diversity of habitats (open and disturbed) colonized by S. oleraceus (Supplementary Table 1). At each site, on each date, we sampled three quadrats (1 m²) along a 20 m linear transect. Quadrats were placed to contain at least one S. oleraceus plant. Within each quadrat, arthropods were collected from plants with a forceps or brush, and stored individually in sterile 2 ml Eppendorf tubes filled with a protective buffer solution. This solution is used to prevent oxidation of polyphenols and polyamines (PCR inhibitors) (see [START_REF] Cruaud | Using insects to detect, monitor and predict the distribution of Xylella fastidiosa : a case study in Corsica[END_REF] for more details). This procedure was repeated for each plant of every plant species present in the quadrat over a period of one hour, to standardize the sampling effort. This period was deemed adapted to represent the biodiversity of the sampled unit, and to allow vagrant insects, potentially disturbed by our arrival, to settle back on their resource plant before sampling. We collected individual specimens except for colonies of aphids, thrips, and egg masses, for which at least five specimens were required to obtain sufficient DNA for analysis. We did not consider pollinators or the soil fauna in this study. Following the collection of each specimen, tools were thoroughly cleaned by successive immersions in 2.5 % bleach solution, water and 96% ethanol, to prevent cross-contamination. At the end of the one-hour insect sampling period, all the plants within the quadrat were collected individually (by cutting the stem at the soil surface), for further dissection.

Back in the laboratory, the plants were identified morphologically, and their organs (stems, leaves, flowers) were dissected to collect endophagous arthropods, which were transferred into tubes as described for the arthropods collected in the field. For each arthropod specimen collected, we identified the plant species from which arthropods were sampled, and recorded the specimen stage and condition (degraded, parasitized), and putative identification (at least taxonomic group, with identification to species level if straightforward). All arthropod samples were frozen at -20°C until DNA analysis. Thus, while plants were identified morphologically, arthropods were identified via molecular technologies. Each plant was transferred to a paper bag and oven-dried at 70°C for 72 h, for the determination of aboveground dry biomass (g) as an estimate of plant abundance per quadrat. Arthropod abundances were determined based on the number of individuals collected per quadrat for each taxon. Sampling was performed for 57 quadrats, over the three sampling dates.

High-throughput DNA metabarcoding

We characterized the interaction network by directly observing plant-arthropod interactions (recording only interactions for which an observation of feeding was verified); while arthropodarthropod interactions were revealed by molecular analysis. We first isolated total DNA from each arthropod individual [START_REF] Cruaud | Using insects to detect, monitor and predict the distribution of Xylella fastidiosa : a case study in Corsica[END_REF]. As presented in Supplementary Figure S1, we then performed metabarcoding on each arthropod sample, with a two-step DNA amplification and high-throughput sequencing method adapted from the procedure described by [START_REF] Galan | Metabarcoding for the parallel identification of several hundred predators and their prey: Application to bat species diet analysis[END_REF]. We sequenced three short COI fragments, with primer combinations and PCR protocols developed elsewhere (HCO forward: [START_REF] Leray | A new versatile primer set targeting a short fragment of the mitochondrial COI region for metabarcoding metazoan diversity: application for characterizing coral reef fish gut contents[END_REF], HCO reverse: Folmer et al. 1994, LEP F. and R.: Brandon-Mong et al. 2015, HEX F. and R.: Marquina et al. 2019), to overcome the problem of the lack of primer universality among arthropods. Error-proof indices for individual sample identification were developed with the high-throughput sequencing process described by [START_REF] Martin | error-proof_indexes v1.0 (Zenodo[END_REF]. The libraries were sequenced with Illumina technology, using a Miseq 2x250 run for date 1 (April), and one lane of Hiseq 3000 each for dates 2 (May) and 3 (June).

The markers for each sample were demultiplexed with CutAdapt v2.3, and all paired-end reads were filtered for minimal length (280 bp), corrected for sequencing errors, and pairs of overlapping reads were merged with the Dada2 v1.12 R package [START_REF] Callahan | DADA2: High resolution sample inference from Illumina amplicon data[END_REF]. A matrix was thus obtained, containing samples as variables and amplicon variant sequences (ASVS) as observations. A variant is a set of identical corrected and merged paired-end reads. We used Qiime2 [START_REF] Bolyen | QIIME 2: Reproducible, interactive, scalable, and extensible microbiome data science[END_REF] with a 2% divergence threshold, to merge ASVS, to decrease their number without the loss of taxonomic information. The summed number of reads for each merged variant for a given arthropod sample was reported as the intersection of samples and ASVS.

Each ASVS was assigned, by BLAST, to a barcoding reference database of cytochrome oxidase subunit I (COI) nucleotide sequences (658 bp) compiled from three different sources and curated by expert analysis. These reference barcodes were retrieved from BOLDSYSTEM [START_REF] Ratnasingham | bold: The Barcode of Life Data System[END_REF], the CBGP -Continental Arthropod collection (Centre de Biologie pour la Gestion des Population, 2019) and a local database specifically designed for this study. Our database contained barcodes of the most frequently encountered species during this sampling campaign (extra-specimens collected) and field surveys (2017-2020) conducted through Europe and North Africa for the search of S. oleraceus natural enemies (see below). In total, these three sources compiled 1 699 995 sequences from 119 299 species available for ASVS assignation. We retained successful assignments to the ranks of species, genus and family, but not those to higher taxonomic levels, because arthropod biology is too variable at higher taxonomic ranks to be informative for our purpose. The assignments obtained for each marker were grouped together in a single table and the numbers of reads were summed by assigned taxon. The resulting file was therefore an interaction matrix in BIOM file format, in which the assigned taxa replaced ASVS. The matrix was curated and manually transformed to obtain an adjacency matrix (in which the observations are sources and the variables are consumers) usable for further network analyses. For each pair of consumer/prey species, occurrence frequencies of interaction were calculated (Supplementary Text 1 and Figure S2).

Assessment of sampling robustness and global network description

We first evaluated the completeness of sampling over the entire sampling campaign, and generated taxon accumulation curves (the 57 quadrats were added in a random order, with 1,000 permutations) for plants and arthropods, using the specaccum function of the R package vegan [START_REF] Oksanen | vegan: Community Ecology Package. R package version 2[END_REF]. We estimated the extrapolated taxonomic richness by calculating the Chao 1 index [START_REF] Chao | Nonparametric Estimation of the Number of Classes in a Population[END_REF] with the specpool function. Likewise, the robustness of sampling for the characterization of interactions was assessed by generating accumulation curves for pairwise interactions. We first generated an accumulation curve including all the types of direct interactions (e.g. plant-herbivores, herbivores-natural enemies, etc.) present in the meta-network (i.e. pooling interactions from all sites). The 57 quadrats were added in a random order, with 1000 permutations. We finally generated a curve focusing on interactions involving S. oleraceus as a source, to evaluate the performance of the sampling design for addressing our objective of establishing an inventory of the arthropods feeding on S. oleraceus, corresponding to candidate BCA. For both curves, we estimated the extrapolated interaction richness with the Chao 1 index [START_REF] Chao | Nonparametric Estimation of the Number of Classes in a Population[END_REF], using the specpool function.

Prior to interaction analyses, a global description of the meta-network (pooling interactions data from all sites) and subnetwork (centred on S. oleraceus related interactions) was performed.

Several metrics were calculated: the number of links (L), the number of nodes (S) (connected and isolated), connectance (C) and link density (LD) [START_REF] Bersier | Quantitative Descriptors of Food-Web Matrices[END_REF][START_REF] Warren | Making connections in food webs[END_REF].

Connectance is the proportion of the possible trophic links actually realized; here cannibalism is not permitted, so C = L/S(S-1). Link density is the mean number of links per taxon, calculated as LD = L/S. We also characterized the taxon assemblage by determining taxonomic richness (i.e. number of taxa) for each trophic level (plants, herbivores and natural enemies).

Selection of candidate biocontrol agents

The selection of candidate BCA was decided according two criteria: a restricted ecological host range and limited interactions with natural enemies. Thus, based on the interactions retrieved from the meta-network, we selected a subnetwork considering only the arthropods having S. oleraceus as a source plant, as well as all their complementary plant resources. We also included natural enemies associated with these herbivores (i.e. parasitoids and predators). We assessed and visualized the specificity of these herbivores, by plotting interactions between herbivores encountered on S. oleraceus and all their complementary resource plants as a grid matrix, in which plants were ordered by their degree of phylogenetic relatedness to S. oleraceus, as defined by the current classification of angiosperms [START_REF] Chase | An update of the Angiosperm Phylogeny Group classification for the orders and families of flowering plants: APG IV[END_REF]. Arthropods were ordered by the increasing generality values (i.e. the number of resources per taxon) characterizing ecological host range. To assess and visualize the dependence of natural enemies on these herbivores, we constructed a second level grid matrix in connexion with the previous, and calculated arthropod vulnerability values (i.e. the number of consumer per taxon). Multipartite network and grid matrices were constructed with igraph R package.

Assessing the contribution of the approach for the biocontrol program To discuss the contribution of the method herein proposed, we used, as a point of reference, a survey performed following classical procedures (sampling, rearing and identification of specimens exclusively collected from S. oleraceus) in the frame of this CBC program (Lesieur et al., in prep). However, we acknowledge that this classical survey covered a longer period of sampling (2017-2020) and a much larger geographical area was prospected (10 countries through Europe and North Africa).

Results

Summary of the molecular results

In total, 2,834 arthropod specimens were collected and analyzed by metabarcoding, to reconstruct the interaction network at a global scale. We obtained DNA sequences and taxonomic assignments for 1,803 of the 2,834 arthropods initially collected (63.6 %). This proportion of exploitable information reached 71% (2,011 specimens) after manual validation of the matrix.

The molecular analysis provided a total of 107,483,410 reads, 19.2% of which were retained after screening with quality filters; we obtained a final dataset of 2,014 COI variant sequences (Supplementary Table 2). Before, manual validation, we observed that a large proportion of the diversity (33% of the families and 40% of the species) was recovered by the use of all markers, the rest being recovered by a combination of two markers, or specifically found with only one marker (Supplementary Figure S3). LEP increased identification rates by 20% for families and 25% for total species, consistent with its widespread use in the research community [START_REF] Brandon-Mong | DNA metabarcoding of insects and allies: an evaluation of primers and pipelines[END_REF]. The other two markers also provided original information, albeit to a lesser extent, at least as far as the number of taxa recovered was concerned, as 15% of the families and 16% of the species would not have been recovered with LEP alone. After data validation, 269 taxa were identified for arthropods, with 84% identified to species level (17 orders, 90 families and 189 genera). While plant taxonomic diversity (relying on morphological identifications) accounted for 132 taxa, 80% of which were classified to species level (25 orders, 29 families and 87 genera) (Supplementary figure S4).

Sampling robustness

The accumulation curve of plants seemed to approach an asymptote, but this was not the case for arthropods (Figure 1). The Chao 1 index indicated an extrapolated taxonomic richness value for plants of 164 taxa (± 12), with 132 taxa actually sampled. By contrast, for arthropods, the extrapolated taxonomic richness value was 442 taxa (± 39), but only 269 taxa were actually sampled. Sampling robustness was high over the entire sampling scheme for plants but sampling efficiency was lower for arthropods. Likewise, we assessed the completeness of pairwise interactions detected over the whole network. We observed a linear increase associated with a Chao1 index of 1245 (± 183) expected interactions, where 350 links were actually reconstructed (Figure 1). However, this is less of an issue for interactions involving S. oleraceus, the focus of the analysis for which this sampling was designed. The accumulation curve in question tended towards an asymptote, with a Chao 1 index of 63 (± 10) expected interactions and 47 interactions sampled. Overall, these results suggest that the sampling effort was adequate for the reconstruction of a unique interaction network maximizing of the proportion of links observed [START_REF] Jordano | Sampling networks of ecological interactions[END_REF].

Meta-network and subnetwork analyses

As presented in Table 1, the complete interaction network (meta-network) consisted of 401 nodes, 241 of which were connected to another node (60%), resulting in 350 links (Supplementary Figure S5). Linkage density and connectance calculated were 1.45 and 0.006, respectively. The meta-network included 60 plants in interaction (46% of the plants collected), 136 herbivores in interaction (74% of the herbivores collected), 35 natural enemies in interaction (49% of the natural enemies collected) comprising 19 parasitoid and 16 predator taxa, and 10 omnivores (feeding at more than one trophic level). The sub-network consisted of 116 nodes and 213 links, and resulting linkage density and connectance were 1.84 and 0.008, respectively (Figure 2). A more detailed description of taxon assemblage composing S. oleraceus subnetwork is provided in the following section.

Identifying candidate biocontrol agents: considering host range and regulation by enemies Analysing S. oleraceus subnetwork, we found 47 herbivorous taxa feeding on the target, including 37 taxa identified to species level. They belonged to five different orders, i.e. Hemiptera (45%), Diptera (25%), Coleoptera (19%), Lepidoptera (0.06%) and Hymenoptera (0.04%), and were distributed in nine different trophic guilds, with the flower bud suckingpiercing guild being the most represented (23%) while the less represented guild corresponded to the chewing guild (2%) (Table 2). Fifteen taxa were collected exclusively from S. oleraceus, and another two taxa were collected from S. oleraceus and Sonchus asper (Figure 3). These taxa are potential BCA (host range apparently restricted to the genus Sonchus, subtribe Sonchinae). Six additional species were detected only on members of the tribe Chicorieae (Aphis craccivora Koch, Ophiomyia cunctata Hendel, Phytomyza lateralis Fallén, Campiglosa producta Loew, L. punctiventris and T. formosa). We identified 38 other plant species as complementary resource plants for the herbivore species collected from S. oleraceus. The generality of these herbivore species ranged from 1 to 18, with Philaenus spumarius L. the most polyphagous of the 47 herbivores species found on S. oleraceus (Figure 3).

The analysis of the subnetwork (Figure 2) also indicated that the herbivores collected on S. oleraceus were a resource for diverse natural enemies. In particular, 19 of the 47 herbivorous taxa collected were attacked by several species of parasitoid (12 species from the family Braconidae, 1 from Figitidae, and 1 from Ichneumonidae) and predators (6 Arachnida species, 1 from Cantharidae, 2 from Coccinnellidae, 3 from Syrphidae and 1 from Orthoptera). Moreover, among the 17 arthropods identified as candidate BCA for their restricted ecological host range, we detected interactions with natural enemies for six of them, one exhibiting interactions with 8 taxa from higher trophic levels (Figure 3).

Eventually, molecular analyses revealed particular patterns of omnivory involving several species from Heteroptera. We distinguished between intermediate omnivores (species feeding on both plants and herbivores, such as members of the Tephritidae and Aphididae), and top omnivores (species feeding on herbivores and natural enemies, such as members of the Syrphidae). The list of the taxa included in the subnetwork and of all the trophic interactions (i.e. the edge list) used to generate Figure 3 are provided in Supplementary Tables 3 and4, respectively.

Assessing the contribution of the approach for the biocontrol program

The analysis of trophic interactions identified 47 taxa feeding on S. oleraceus, 37 of which were identified to species level. Nineteen of these 37 species had already been sampled in classical field surveys, the other 18 species being newly reported as herbivores of S. oleraceus in this CBC program (Table 2).

Discussion

The combination of observation and molecular data performs well for the characterization of interactions

The metabarcoding approach used made it possible to target a broad diversity of taxa (90 arthropod families, with identification to species level of 84% of the variants), as expected with the use of multiple markers (Alberdi et al., 2018;[START_REF] Creedy | Toward accurate species-level metabarcoding of arthropod communities from the tropical forest canopy[END_REF][START_REF] Marquina | New mitochondrial primers for metabarcoding of insects, designed and evaluated using in silico methods[END_REF]. The combination of taxonomic assignments with subsequent observational data and information available in the literature was essential: 1) to validate the trophic links (predation and parasitism) and 2) to complement the identification in cases of failed amplification or taxonomic assignment, as advocated in other contexts [START_REF] Derocles | Chapter One -Biomonitoring for the 21st Century: Integrating Next-Generation Sequencing Into Ecological Network Analysis[END_REF][START_REF] Wirta | Complementary molecular information changes our perception of food web structure[END_REF]. In the meta-network, 60% of the taxa interacted, suggesting that our methods performed very well for the reconstruction of interactions. More specifically, interaction detection rates obtained for plantherbivores and herbivores-natural showed higher values than those usually reported in comparable contexts [START_REF] Braukmann | Testing the Efficacy of DNA Barcodes for Identifying the Vascular Plants of Canada[END_REF][START_REF] Clare | Molecular detection of trophic interactions: emerging trends, distinct advantages, significant considerations and conservation applications[END_REF][START_REF] Erickson | Reconstructing a herbivore's diet using a novel rbcL DNA mini-barcode for plants[END_REF][START_REF] García-Robledo | Tropical plant-herbivore networks: reconstructing species interactions using DNA barcodes[END_REF][START_REF] Roslin | The use of DNA barcodes in food web construction-terrestrial and aquatic ecologists unite![END_REF]. The high rate of interaction reported here for herbivorous arthropods can be explained by our decision to focus on intensive plant dissection and morphological determination. Retaining feeding interactions only after verification reduced the risk of false positives, over-estimating species interactions, related to the use of co-occurrence data (i.e. tourist insects on plants rather than actual trophic links) [START_REF] Zhu | Seeing is believing? Comparing plant-herbivore networks constructed by field co-occurrence and DNA barcoding methods for gaining insights into network structures[END_REF]. For arthropods, the rather low rate of natural enemies positive for preys can be multifactorial; e.g. mismatch between the primer pairs used and the prey species, low sequencing depth given the DNA yield ratio between consumer and prey, and degradation of DNA from consumed preys [START_REF] Hosseini | Factors affecting detectability of prey DNA in the gut contents of invertebrate predators: a polymerase chain reaction-based method[END_REF][START_REF] Macías-Hernández | Molecular gut content analysis of different spider body parts[END_REF][START_REF] Sheppard | Advances in molecular ecology: tracking trophic links through predator-prey food-webs[END_REF][START_REF] Sheppard | Infiltration by alien predators into invertebrate food webs in Hawaii: a molecular approach[END_REF].

A complementary inventory of herbivores feeding on S. oleraceus

The 47 taxa collected from S. oleraceus, covered a wide range of trophic guilds. This evidenced that the method herein employed did not bias the selection towards a particular trophic guild but allows the detection of herbivores exhibited diversified feeding habits. The selection of one or several BCA from these trophic guilds could offer a good complementarity of actions [START_REF] Buccellato | Post-release evaluation of a combination of biocontrol agents on Crofton weed: testing extrapolation of greenhouse results to field conditions[END_REF]. Although the guild of flower head sucker-piercer was the richest, further dedicated experiments would be necessary to assess whether those candidate agents actually provide the best regulation action of S. oleraceus [START_REF] Morin | Review of approaches to evaluate the effectiveness of weed biological control agents[END_REF].

Several elements indicated that our approach is a good complement of the classical survey. For example, we provided a more detailed description of pollen chewers (only Brassicogethes aeneus (Fabricius, 1775) was recorded in the classical survey). For Diptera, the Tephritidae flies, Tephritis cometa Loew and T. vespertina Loew, are newly recorded. A phylogeny of this taxa based on the CO1 barcode showed that those species are closely related to T. formosa Loew [START_REF] Smit | Half of the European fruit fly species barcoded (Diptera, Tephritidae); a feasibility test for molecular identification[END_REF], and we cannot, therefore, rule out possible molecular misidentification due to the short CO1 barcode used or host race differentiation, as frequently observed in this group [START_REF] Diegisser | Phylogeographic patterns of host-race evolution in Tephritis conura (Diptera: Tephritidae)[END_REF]. During the classical survey, despite intensive collections, the Tephritidae species Campiglossa producta has been sampled only on S. oleraceus from the Canary Islands.

The intensive sampling efforts put in the present study on S. oleraceus led to the collection of rare C. producta specimens (only 15 specimens from two sites) in the continental bioclimatic region in France. This highlights that, despite a reduced area prospected and a limited time frame, some rare species were sampled and their ecological host range described. We acknowledge we missed some species occurring later in the season or out of the sampled area. For example, Cystiphora sonchi Vallot (Diptera: Cecidomyiidae) was collected in the classical surveys and passed specificity tests [START_REF] Lesieur | Making host specificity testing more efficient: Exploring the use of abridged test plant lists[END_REF], but was not sampled from S. oleraceus in the three bioclimatic regions from April to June. This lack of detection in our sampling campaign was expected, as rates of infestation with this species peak in summer [START_REF] Rizzo | On two Italian gall midges (Diptera Cecidomyiidae) and their parasitoids[END_REF].

Hence, this study should be regarded as a complement of usual procedures. However, with the rapid development of molecular technologies and associated drop in price [START_REF] Kennedy | High-throughput sequencing for community analysis: the promise of DNA barcoding to uncover diversity, relatedness, abundances and interactions in spider communities[END_REF], we believe that this approach will be soon applicable at larger sampling scales.

Selection of BCA First criterion: a restricted ecological host range

Based on the present results, 15 of the 47 herbivores feeding on S. oleraceus seemed to have an ecological host range restricted to S. oleraceus, and another two taxa appeared to be restricted to the genus Sonchus. All these taxa recovered from the ecological network are, thus, of particular interest as candidate BCA. However, contradictions were observed between the ecological host range described by network analysis and published findings or specificity test results (Lesieur et al., in prep). These discrepancies may be due to insufficient sampling for the recording of species interactions (as shown by accumulation curve on total interactions) or to the presence of cryptic host races or cryptic species that have yet to be deciphered. Additional studies will be required to characterize the ecological host ranges of these species further. In particular, two species, Liriomyza sonchi Hendel and Ensina sonchi L., were found associated with S. oleraceus and S. asper and were, therefore, considered to be candidate BCA because these plants are both invasive weeds in Australia [START_REF] Cullen | Biological Control of Weeds in Australia[END_REF]. However, a wider range of food resources has been reported in literature for these two species (Table 2). Conversely, the promising galling insect T. formosa passed specificity tests during the traditional phase of the CBC program. It was found to be restricted to the genus Sonchus, contrary to the results reported here, as we found T. formosa on Crepis vesicaria L. (belonging to Cripidinae, the same tribe but a different subtribe to S. oleraceus.). This plant was not indented to be tested as a potential food plant for T. formosa, and these results therefore highlight the complementarity of the ecological network approach for clarifying herbivore host range.

Moreover, in the interaction network, C. producta was identified on three different plant species from the Chicorieae tribe. This species, found only on S. oleraceus in the Canary Islands during the classical surveys, was considered a promising BCA for testing. The results presented here indicate that its host range would not be compatible with its use as a BCA, potentially leading to its exclusion from the list of candidate BCA. This example shows how the network developed here is complementary to classical procedures, making it possible to narrow down the list of candidate BCA to be tested.

The same applies to Cheilosia latifrons Zetterstedt, a species collected in the classical survey. In our study, we did not sample this species on S. oleraceus, but the meta-network indicated it was collected from S. asper and Picris echioides (L.), revealing its oligophagous dietary behavior. However, little is known about the biology and the host plants of C. latifrons [START_REF] Schmid | Eiablage von Cheilosia latifrons (Zetterstedt, 1843)(Díptera, Syrphidae) an Leontodón autumnalis L. Volucella[END_REF] and its taxonomy seems to be unsolved, calling into question the existence of a species complex defined on the basis of host plant use [START_REF] Speight | Species accounts of European Syrphidae (Diptera)[END_REF]. We also observed discrepancies for specimens from Cynipidae that appeared to be generalist herbivores in the interaction network (associated with both S. oleraceus and Carduus pycnocephalus Spreng.), whereas subsequent analysis of the variants assigned to Cynipidae indicated a genetic structure more consistent with multiple cryptic species potentially specializing on the host plants from which they were collected. One species from Cynipidae is a known stem galler of S. oleraceus: Aulacidea follioti Barbotin (Bladmineerders Online database, 2020). However, this species is not yet present in any of the barcoding databases used here and could therefore only be assigned to family level. Further prospections to collect other Cynipidae specimens and rear them to adulthood would be required to confirm this identification.

Second criterion: limited interactions with natural enemies By using metabarcoding to reconstruct interactions between arthropods, we were able to detect a wide range of parasitoids from their herbivore hosts, and some predators. We detected omnivorous dietary behavior in several groups from Heteroptera. Opportunistic predation through carnivory is common in Lygaeidae [START_REF] Burdfield-Steel | The evolutionary ecology of the Lygaeidae[END_REF] to supplement the low levels of protein supplied by plants. Carnivory has also been reported in Miridae [START_REF] Wheeler | Miridae): Pests, Predators, Opportunists[END_REF] and sometimes leads to intraguild predation interactions. We found that both Syrphidae (Diptera) and Miridae (Heteroptera) fed on aphid species, but we also revealed that mirids could prey upon syrphids, as already demonstrated in arena experiments [START_REF] Fréchette | Intraguild predation between syrphids and mirids: who is the prey? Who is the predator?[END_REF]. Members of the Lygaeidae and Miridae were also found to be able to access and feed on larval stages of several Tephritidae species whilst inside the flower heads of S. oleraceus. This interaction does not seem to have been observed before and provides insight useful not only for the CBC program against S.

oleraceus, but also with direct implications for other biological control programs, particularly those involving the conservation biological control of insect pests.

More specifically, among candidate BCA exhibiting a restricted ecological host range, some were associated to an important diversity of natural enemies, and should be considered of lower priority for testing (i.e. Ensina sonchi). We suggest that priority should be given to agent predated by a minimum parasitoid and predator species, and, when they exist, to an agent that possesses enemies from the most distant taxonomical group from those occurring in the range of introduction [START_REF] Ollivier | Characterizing ecological interaction networks to support risk assessment in classical biological control of weeds[END_REF]. It has been shown that newly created interactions between hosts and parasitoids in the introduced range are predictable based on the realised interactions in native range [START_REF] Paynter | Making weed biological control predictable, safer and more effective: perspectives from New Zealand[END_REF][START_REF] Veldtman | Creating novel food webs on introduced Australian acacias: indirect effects of galling biological control agents[END_REF]. Further steps in this program would consist in investigating the diversity of natural enemies occurring in the range of invasion to anticipate new potential interactions and refine BCA choice.

Conclusion

We demonstrate here the potential of network ecology for characterizing candidate BCA and their ecological interactions in the field. This characterization clearly benefited from the use of complementary approaches (morphological and molecular analyses) to identify plant/arthropod and arthropod/arthropod interactions and provided a solid framework for the establishment of an inventory of herbivores feeding on the target weed, their realized host range and interactions with natural enemies. Avenues for further investigation have been identified and in-depth studies are now required. The strength of this approach also lies in its capacity to screen field host ranges for multiple herbivore species simultaneously, without the need for as many tests as species. This potential to narrow down the list of candidate BCA for testing should help to save both time and money. Finally, in addition to the potential value of ecological network analysis to the CBC targeting Common Sowthistle in Australia, the data reported here are potentially useful for other future programs.

Authors' contributions

Conceptualization: M. Ollivier, V. Lesieur, M. S. Tixier, J.-F. Martin Methodology: M. Ollivier, J. Tavoillot, F. Bénetière, V. Lesieur, J.-F. Martin Tables Table 1: Global description of the meta-network and subnetwork centred on S. oleraceus. Metrics measured are the total number of nodes (in brackets is the number of nodes connected in the network), the number of links, the linkage density and the connectance of the network. For each network is also presented the number of species per trophic level (in brackets is the number of species in interaction with another species). Omnivorous species are regarded as natural enemies. 3 and4, respectively.

Meta-network

Supplementary Texts

Text 1: Curation and transformation of the interaction matrix

In order to produce an interaction matrix, in which observations are resources and variables are consumers, the original dataset should be checked. First, for each sample, the taxa having the highest number of reads were set as the consumer. This assumption was confirmed in most of the cases (the specimen collected in the tube was the consumer and produced the highest number of reads). However due to the simultaneous detection of possible preys and parasitoids from the specimen sampled and the contrasted amplification rates between taxonomic groups (Alberdi et al. 2018), this was not always the case. Consequently, a second step of data manual validation was required on each sample (Supplementary Figure S2) in order to 1) verify whether the consumer was correctly defined according to the observation from the field, 2) add the source plant for herbivorous arthropods, 3) define the status of the taxa simultaneously amplified with the consumer (prey, parasitoid or contamination). The next step consisted in removing cannibalistic interactions (i. e. delete reads at the intersection of a same source and consumer species). Indeed, DNA-based methods are not suitable to disentangle predator and prey sequences [START_REF] Taberlet | Environmental DNA: For Biodiversity Research and Monitoring[END_REF]. As metabarcoding neither provides quantitative information about the trophic links [START_REF] Deagle | Counting with DNA in metabarcoding studies: How should we convert sequence reads to dietary data?[END_REF], the numbers of reads were converted to binary information, then occurrence frequencies of interactions between a pair of consumer/prey species were calculated as the percentage of consumer samples positive for the prey over the total number of consumer sampled in the dataset. All analyses were carried out under R version 3.5.2 (R Core Team 2018).

Figure 1 :

 1 Figure 1: Accumulation curves representing a) species richness in plants and arthropods and b) pairwise interactions from the meta-network and focusing on interactions involving S. oleraceus. Curves were constructed with 1000 random resampling events over the 57 quadrats analyzed along the sampling campaign (Spring 2018).

Figure 2 :

 2 Figure 2: Multitrophic network reconstructed from S. oleraceus (dark green node) and 38 other plants (light green nodes) used by S. oleraceus herbivores (red nodes represent putative specialist herbivores, and pink nodes correspond to herbivores feeding on other species in addition to S. oleraceus). Plants are ordered by their phylogenetic relatedness to S. oleraceus. Natural enemies of herbivore species are represented by dark blue nodes. Nodes at intermediate levels (beige and light blue) correspond to omnivorous species identified by molecular analyses. The width of edges reflects frequencies of interactions between pairs of species and edges coloured in red emphasize interactions involving potential candidate BCA. The list of taxa corresponding to each node and the edge list are provided in Supplementary Tables3 and 4, respectively.

Figure 3 :

 3 Figure 3: Interaction matrix between herbivores sampled on Sonchus oleraceus and their resource plants, indicating the ecological host range of the herbivores, as defined by intense field sampling in France, during spring 2018. Plants are ordered by their phylogenetic relatedness to S. oleraceus and arthropods are ordered by increasing generality values (i.e. the number of resources per species).Red rectangle highlights the 17 species considered as candidate BCA for their restricted eclogical host range. The second level matrix diplays interactions between candidate BCA and their natural enemies. Interaction are represented as semi-quantitative information, via occurences frequencies of interactions culculated for each species pair.

Figure S1 :

 S1 Figure S1: Overview of the entire process for the reconstruction of interactions between plants, herbivores and natural enemies.

Figure S3 :

 S3 Figure S3: Venn diagram of taxonomic coverage obtained with the three COI markers used in this study, represented through a) arthropod family diversity and b) arthropod species diversity.

Table 2 :

 2 Herbivores from Sonchus oleraceus (SO) collected and identified through intense field sampling in three bioclimatic regions (Semi-oceanic, Mediterranean and Continental) in France in spring 2018. The field host range of these herbivores was defined by network analysis (S -SO: Specific to SO, S -S: Specific to the genus Sonchus, S -C: Specific to the tribe Cichorieae, G: Generalist, ?: Unknown host range). The column 3-years surveys refers to the species collected following classical procedures of biocontrol program (see Materials and Methods section).

	Order / Family	Species	Trophic guild	3-years survey	Field host range	Literature host range (Reference)
	Coleoptera
	Cerambicydae Agapanthia cardui (Linnaeus, 1767)	Stem boring	No	G
	Cerambicydae Agapanthia pannonica Kratochvil, 1985	Stem boring	No	S -SO	Polyphagous (Wang, 2017)
	Curculionidae sp.	Crown/stem boring	_	G
	Curculionidae Lixus punctiventris (Boheman, 1835)	Stem boring	Yes	S -C
	Dasytidae	Dasytes tristiculus Mulsant & Rey, 1868	Flower/pollen chewing	No	S -SO	Polyphagous (Diputacio Barcelona, 2019)
	Dasytidae	Psilothrix viridicoeruleus (Geoffroy, 1785)	Flower/pollen chewing	No	G
	Oedemeridae	Oedemera crassipes Ganglbauer, 1881)	Flower/pollen chewing	No	S -SO	Polyphagous (Vázquez Albalate, 2002)
	Oedemeridae	Oedemera flavipes (Fabricius, 1792)	Flower/pollen chewing	No	G
	Nitidulidae	Brassicogethes aeneus (Fabricius, 1775)	Flower/pollen chewing	Yes	G
	Diptera
	Agromyzidae	sp.	Leaf mining	_	G
	Agromyzidae	Liriomyza sonchi Hendel, 1931	Leaf mining	Yes	S -S	Oligophagous (Benavent-Corai, 2005)
	Agromyzidae	Ophiomyia cunctata (Hendel, 1920)	Leaf mining	Yes	S -C
	Agromyzidae	Phytomyza horticola Goureau, 1851	Leaf mining	Yes	G
	Agromyzidae	Phytomyza lateralis Fallén, 1823	Leaf mining	Yes	S -C
	Tephritidae	sp.	Flower bud galling/seed feeding _	S -SO	Oligophagous (White, 1988)
	Tephritidae	Campiglossa producta (Loew, 1844)	Flower bud galling/seed feeding Yes	S -C
	Tephritidae	Ensina sonchi (Linnaeus, 1767)	Flower bud galling/seed feeding Yes	S -S	Oligophagous (White, 1988)
	Tephritidae	Tephritis sp.	Flower bud galling/seed feeding _	G
	Tephritidae	Tephritis cometa (Loew, 1840)	Flower bud galling/seed feeding No	S -SO	Oligophagous (Bladmineerders Online database, 2020)
	Tephritidae	Tephritis formosa (Loew, 1844)	Flower bud galling/seed feeding Yes	S -C

Table 4 :

 4 List of interactions included in the sub-network(Figure 2)

	Interaction Resource	Consumer	Occurrence frequency
	1	Sonchus asper	Agromyzidae	20
	2	Sonchus asper	Campiglossa producta	13
	3	Sonchus asper	Ensina sonchi	9
	4	Sonchus asper	Hyperomyzus lactucae	17
	5	Sonchus asper	Liorhyssus hyalinus	34
	6	Sonchus asper	Liriomyza sonchi	14
	7	Sonchus asper	Lixus punctiventris	2
	8	Sonchus asper	Ophiomyia cunctata (beckeri)	8
	9	Sonchus asper	Philaenus spumarius	9
	10	Sonchus asper	Phytomyza horticola (syngenesiae)	15
	11	Sonchus asper	Psilothrix viridicoeruleus	9
	12	Sonchus asper	Tephritis formosa	7
	13	Sonchus asper	Tephritis sp	11
	14	Sonchus oleraceus	Agapanthia cardui	15
	15	Sonchus oleraceus	Agapanthia pannonica	50
	16	Sonchus oleraceus	Agromyzidae	52
	17	Sonchus oleraceus	Aphididae	100
	18	Sonchus oleraceus	Aphis craccivora	50
	19	Sonchus oleraceus	Aphis fabae	71
	20	Sonchus oleraceus	Brassicogethes aeneus	21
	21	Sonchus oleraceus	Campiglossa producta	60
	22	Sonchus oleraceus	Cicadella viridis	100
	23	Sonchus oleraceus	Closterotomus norvegicus	100
	24	Sonchus oleraceus	Cnephasia stephensiana	33
	25	Sonchus oleraceus	Curculionidae	12
	26	Sonchus oleraceus	Cynipidae	50
	27	Sonchus oleraceus	Dasytes tristiculus	8
	28	Sonchus oleraceus	Dolycoris baccarum	50
	29	Sonchus oleraceus	Ensina sonchi	90
	30	Sonchus oleraceus	Hecatera dysodea	100
	31	Sonchus oleraceus	Hyalopterus pruni	100
	32	Sonchus oleraceus	Hyperomyzus lactucae	69
	33	Sonchus oleraceus	Lepidargyrus ancorifer	22
	34	Sonchus oleraceus	Liorhyssus hyalinus	50
	35	Sonchus oleraceus	Liriomyza sonchi	79
	36	Sonchus oleraceus	Lixus punctiventris	96
	37	Sonchus oleraceus	Lygaeus equestris	11
	38	Sonchus oleraceus	Macrosiphum rosae	100
	39	Sonchus oleraceus	Miridae	38

Acknowledgments

We are particularly grateful to A Coeur d'acier, G Delvare, B Derepas, J Harran, B Michel, A Migeon, E Pierre, JM Ramel and JC Streito for providing their taxonomic expertise for arthropod identification for construction of the local database. We also thank G Fried for assistance with the identification of some plant specimens. We also warmly thank all our coworkers involved in field sampling: P Audiot, M Corbin, R Guilhot, A Loiseau, L Olazcuaga, T Thomann, and N Vieira. This project is supported by funding from the Australian Government Department of Agriculture and Water Resources, as part of its Rural R&D for Profit programme, through AgriFutures Australia (Rural Industries Research and Development Corporation) (PRJ--010527).

Data acquisition: M. Ollivier, J. Tavoillot, F. Bénetière, V. Lesieur, J.-F. Martin Data analysis: M. Ollivier Supervision: M. S. Tixier, J.-F. Martin, Writing original draft: M. Ollivier Writing review and editing: M. Ollivier, J. Tavoillot, V. Lesieur, M. S. Tixier, J.-F. Martin

Supplementary materials

Supplementary Texts Supplementary Tables Supplementary Figures References Supplementary Tables Table 1: Locations of the different sites prospected per sampling session in the three bioclimatic regions.


Aymen Kammoun
email: kammoun@i3s.unice.frwww.i3s.unice.fr

Frédéric Payan
email: fpayan@i3s.unice.frwww.i3s.unice.fr

Marc Antonini

Bit Allocation For Spatio-temporal Wavelet Coding Of Animated Semi-regular Meshes

Keywords: Spatio-temporal decomposition, Wavelet filtering, Animation, Semi-regular meshes, Bit Allocation, Coding

In this paper, we propose a compression scheme for animated semi-regular meshes. This scheme includes a spatio-temporal wavelet filtering to exploit the coherence both in time and space. In order to optimize the quantization of both spatial and temporal wavelet coefficients, the proposed compression scheme also includes a model-based bit allocation. The experimental results show that this approach significantly improves the compression performances, when comparing with previous similar approaches.

Introduction

Typically, an animation is represented by a sequence of triangular meshes sharing the same connectivity at any frame. Even if the static triangulation has the advantage to strongly reduce the amount of information needed to represent such data, compression is useful when the number of vertices increases, because of the amount of geometric information.

During the last decade, most of works focused on the compression for sequences of irregular meshes. When observing the state of the art for such data, we can note techniques using spatial and/or temporal prediction [START_REF] Ibarria | Dynapack: space-time compression of the 3D animations of triangle meshes with fixed connectivity[END_REF][START_REF] Yang | Compression of 3D triangle mesh sequences based on vertex-wise motion vector prediction[END_REF]; principal component analysis (PCA) [START_REF] Alexa | Representing animations by principal components[END_REF][START_REF] Karni | Compression of soft-body animation sequences[END_REF][START_REF] Vasa | Coddyac: Connectivity driven dynamic mesh compression[END_REF]; motion-based segmentation [START_REF] Boulfani | Motion-based mesh clustering for mcdwt compression of 3d animated meshes[END_REF][START_REF] Mamou | Two optimizations of the mpeg-4 famc standard for enhanced compression of animated 3d meshes[END_REF]; waveletbased analysis tool [START_REF] Lopes | Wavelet compression and transmission of deformable surfaces over networks[END_REF][START_REF] Guskov | Wavelet compression of parametrically coherent mesh sequences[END_REF][START_REF] Payan | Temporal wavelet-based geometry coder for 3d animations[END_REF][START_REF] Cho | 3d dynamic mesh compression using wavelet-based multiresolution analysis[END_REF].

Among the wavelet-based methods, the pioneering coder for sequences of irregular meshes was based only on temporal wavelet filtering [START_REF] Lopes | Wavelet compression and transmission of deformable surfaces over networks[END_REF]. The reason is that the geometry of an animation has a regular temporal sampling (same connectivity at any frame), while having an irregular "spatial" sampling (the frames are irregular meshes). And it is well-known that wavelets are more efficient on regular grids than irregular ones (in image coding for instance). Recently, another coder using temporal wavelet filtering but implemented in lifting scheme was proposed [START_REF] Payan | Temporal wavelet-based geometry coder for 3d animations[END_REF]. In parallel, a motion-compensated temporal wavelet filtering has been also developed in [START_REF] Boulfani | Motion-based mesh clustering for mcdwt compression of 3d animated meshes[END_REF]. The main idea of all these coders is to apply a one-dimensional wavelet filtering on the successive positions of each vertex.

Despite the challenge of applying wavelets on irregular grids, some coders using irregular spatial wavelet was proposed. For instance, Guskov and Khodakovsky have developed in [START_REF] Guskov | Wavelet compression of parametrically coherent mesh sequences[END_REF] a coder for parametrically coherent mesh sequences, using an irregular wavelet-based decomposition. Unfortunately, from an analysis point-of-view, all these methods tend to be suboptimal since wavelets are not fully exploited both in spatial and temporal dimension. With this outlook, Cho et al. recently proposed in [START_REF] Cho | 3d dynamic mesh compression using wavelet-based multiresolution analysis[END_REF] to expand the wavelet-based coder for static irregular meshes of [START_REF] Valette | A wavelet-based progressive compression scheme for triangle meshes : Wavemesh[END_REF] to the time-varying surfaces.

In geometry compression, it has been shown that combining remesher and wavelets based on semi-regular sampling leads to efficient techniques for compressing static shapes [START_REF] Khodakovsky | Progressive geometry compression[END_REF][START_REF] Khodakovsky | Normal mesh compression[END_REF][START_REF] Payan | An efficient bit allocation for compressing normal meshes with an error-driven quantization[END_REF]. The semi-regular meshes have the advantage to present a multiresolution structure which facilitates wavelet filtering (to exploit the spatial coherence), and also visualization, levels of details, progressive transmission and so on. In order to have the same properties for animations, similar approaches have been lately proposed for animations defined by sequences of meshes with static or independent connectivity [START_REF] Payan | Remeshing and spatio-temporal wavelet filtering for 3d animations[END_REF][START_REF] Yang | Semi-regular representation and progressive compression of 3d dynamic mesh sequences[END_REF]. The resulting semiregular animations finally have the same relevant properties as the static meshes and also a temporal coherence which can be exploited to perform efficient compression.

We propose in this paper a spatio-temporal wavelet-based coding scheme for sequences of semi-regular meshes sharing the same connectivity. The proposed coder includes an improved version of the bit allocation process presented in [START_REF] Payan | Remeshing and spatio-temporal wavelet filtering for 3d animations[END_REF], which optimizes the quantization of the temporal wavelet coefficients but also the spatial ones (contrary to [START_REF] Payan | Remeshing and spatio-temporal wavelet filtering for 3d animations[END_REF]). The key contribution of this paper is a formulation of the distortion criterion which takes into account the quantization errors of all the data (temporal wavelet coefficients, spatial wavelet coefficients and low frequency signal).

The rest of this paper is organized as follows. Section 2 gives an overview of the proposed method. Section 3 shows the principle of a spatio-temporal wavelet filtering. Section 4 shows the bit allocation for a spatio-temporal decomposition. Section 5 shows some results for several animations, and we conclude in section 6.

Overview of the proposed approach

Let us consider an animated semi-regular mesh F, represented by a sequence of T meshes sharing the same connectivity: F = {f 1 , f 2 , ..., f t , ...f T }. Each mesh f i corresponds to one frame of the animation. Semi-regular animations can be obtained by using a remeshing technique such as [START_REF] Payan | Remeshing and spatio-temporal wavelet filtering for 3d animations[END_REF][START_REF] Yang | Semi-regular representation and progressive compression of 3d dynamic mesh sequences[END_REF]. Here is the outline of our method, illustrated by figure 1. -Spatio-temporal Analysis: The animation F has a regular sampling in time but also A semi-regular one in "space". So, a temporal wavelet-based filtering is first applied, followed by the spatial butterfly-based lifting scheme [START_REF] Khodakovsky | Progressive geometry compression[END_REF]. Hence, we obtain a spatio-temporal decomposition of the input animation (see figure 2): several sets of temporal and spatio-temporal highfrequency details (or wavelet coefficients), and a short sequence of coarse meshes (low frequency signal ). The list of triangles of one sole coarse mesh represents the whole connectivity information.

-Bit Allocation: this process optimizes the rate-distortion trade-off relative to the data quantization. The objective is to efficiently dispatch the bits across the subbands of wavelet coefficients, according to their influence on the quality of the reconstructed mesh sequence. With this outlook, this process will compute the optimal quantization step for each subband, for one usergiven total bitrate.

-Entropy coding: once quantized, the coefficients are entropy coded with the context-based arithmetic coder of [START_REF] Payan | Temporal wavelet-based geometry coder for 3d animations[END_REF].

-Connectivity coder: to reconstruct the compressed data after storage and/or transmission, the connectivity information must be preserved. For this purpose, the list of triangles is encoded in parallel with [START_REF] Touma | Triangle mesh compression[END_REF] and included in the compressed bitstream.

Spatio-temporal decomposition

To create a spatio-temporal decomposition of a semi-regular animation, we choose to apply successively a temporal wavelet filtering and a spatial one. These two filterings are based on lifting scheme implementations [START_REF] Sweldens | The lifting scheme: A custom-design construction of biorthogonal wavelets[END_REF].

Temporal wavelet filtering

The principle of our temporal wavelet filtering is to consider each vertex trajectory as a one-dimensional (1D) signal [START_REF] Payan | Temporal wavelet-based geometry coder for 3d animations[END_REF]. Therefore, a classical 1D lifting scheme is applied on the successive positions of each vertex in parallel. After this step, we obtain two subsets: a high frequency subband H t 1 defined by T /2 detail frames (temporal wavelet coefficients):

H t 1 = {h t 1,0 , h t 1,1 , ...,

Spatial wavelet filtering

After applying the temporal wavelet filtering, we apply the spatial wavelet filtering based on the butterfly lifting scheme [START_REF] Khodakovsky | Progressive geometry compression[END_REF] on all the frames of the low frequency sequence L t Nt . For each frame l t Nt,k of this sequence, this step gives N s subbands of spatial wavelet coefficients h s l t r s ,k (where r s represents the spatial resolution level) and a coarse version l s l t Ns,k of the given frame. The principle of this decomposition is illustrated by figure 2. In this figure, we apply 2 levels of temporal decomposition, and 2 levels of spatial decomposition on the animated Dolphin.

… … … … … … t h 0 , 1 t h 1 , 1 t h 2 , 1 t h 0 , 2 t hl 0 , 1 t hl 1 , 1 t hl 0 , 2 t hl 1 , 2
Fig. 2. Spatio-temporal decomposition based on wavelet filtering: 2 levels of temporal decomposition giving the temporal wavelet coefficients (in red and orange), and 2 levels of spatial decomposition giving the spatial ones (in blue). The meshes in the dotted rounded box (in green) represent the coarse sequence (spatio-temporal low frequency signal).

Given the resulting output multiresolution structure of this analysis stage, we now have to deal with the necessary bit allocation process which will optimize the quantization of the different subbands.

Bit allocation

Given a multiresolution structure, we can not use the same quantizer for all the subbands since the energy is not uniformly dispatched between them. We recall that once transformed, the relevant data (from a coding point of view) is most of times concentrated in the lower frequency subbands, while having the fine and perhaps negligible details in the highest frequency subbands. Therefore we generally use a bit allocation process, to compute the best quantizer for each subband, in order to obtain the best trade-off between rate (compressed file size) and quality of the reconstructed data (distortion).

In our case, the general purpose of the bit allocation process is precisely to determine the best set of quantization steps {q s i,j , q t i,j } -where q s i,j and q t i,j are respectively the quantization steps for the spatial subbands and for the temporal subbands-that minimizes the reconstruction error D T ({q s i,j , q t i,j }) for a given rate R target . This can be formulated by the following problem:

(P) minimize D T {q s i,j , q t i,j } with constraint R T {q s i,j , q t i,j } = R target . (1)
Before solving this problem, we first must formulate the reconstruction error D T .

Note: To know how solving such an optimization problem and have more details, please refer to [START_REF] Payan | An efficient bit allocation for compressing normal meshes with an error-driven quantization[END_REF].

Distortion criterion

Temporal decomposition Let us consider a temporal wavelet coder, using one level of decomposition (see figure 3). In this case, the reconstruction error 3. Model of a temporal wavelet coder used to link the quantization effects with the reconstruction error after decoding. (h0,1, g0,1) and (h1,1, g1,1) are respectively the lowpass analysis/synthesis filters and the high-pass ones. Q represents the quantization process. ↓ 2 and ↑ 2 are respectively the downsampling and upsampling stages. D T can be written as the weighted sum of the MSE (σ t 01 ) 2 and (σ t 11 ) 2 relative to the quantization of the low frequency signal and of the wavelet coefficients [START_REF] Usevitch | Optimal bit allocation for biorthogonal wavelet coding[END_REF]:

+ S S 1 s 0 s 1 ŝ0 ŝ 1 , 0 g 1 , 1 g 1 , 0 h 1 , 1 h 2 2 2 2 Q Q Fig.
D T = 1 2 w t 01 (σ t 01 ) 2 + 1 2 w t 11 (σ t 11 ) 2 , ( 2
)
where w t ij are the weights due to the non-orthogonality of the temporal wavelet filters (i = 0 gives a low frequency subband, i = 1 a high frequency subband and j is the resolution level). Those weights are given by:

w t 01 = 2 N tr (G t 01 ) T G t 01 (3)
and

w t 11 = 2 N tr (G t 11 ) T G t 11 , ( 4
)
where N is the number of vertices and G t ij is the matrix of reconstruction relative to the synthesis filter g t ij . Let us consider now a temporal wavelet coder with N t levels of decomposition. By using (2), we find that the reconstruction error at any level of decomposition (n -1) (2 ≤ n ≤ N t ) can be given by:

(σ t 0(n-1) ) 2 = 1 2 w t 0n (σ t 0n ) 2 + 1 2 w t 1n (σ t 1n ) 2 . ( 5
)
Since we use the same synthesis filter for each level of decomposition,

g 0i = g 01 and g 1i = g 11 ∀i ∈ [2, N t ], (6)
and consequently, the weights are also similar whatever the level:

w t 0i = w t 01 = w t lf and w t 1i = w 11 = w t hf ∀i ∈ [2, N t ]. ( 7
)
By combining the new weights of ( 7), ( 2) and ( 5), we finally obtain:

D T = 1 2 Nt (w t lf ) Nt (σ t 0Nt ) 2 + j=Nt j=1 1 2 j (w t lf ) j-1 w t hf (σ t 1j ) 2 . ( 8
)

Spatial decomposition

The distortion for N s levels of spatial decomposition of a given triangular mesh can be estimated by [START_REF] Payan | An efficient bit allocation for compressing normal meshes with an error-driven quantization[END_REF]:

(σ S ) 2 = (w s lf ) Ns (σ s 0Ns ) 2 + j=Ns j=1 (w s lf ) j-1 w s hf (σ s 1j ) 2 , ( 9
)
where w s lf and w s hf are the weights due to the non-orthogonality of the spatial wavelet transform for the low and high frequency subbands, respectively.

Spatio-temporal decomposition Now, we can generalize the distortion criterion for a spatio-temporal decomposition with N t levels of temporal decomposition and N s levels of spatial decomposition. For each frame of the temporal low frequency subband, we obtain N s spatial subbands. So, the number of spatial subbands is N BF N s , where N BF is the number of frames in the temporal low frequency subband. By replacing (σ t 0Nt ) 2 by ( 9) in ( 8), we finally obtain:

D T = 1 2 N t (w t lf ) Nt f =NBF f =1 1 NBF (w s lf ) Ns (σ s 0Ns ) 2 + j=Ns j=1 (w s lf ) j-1 w s 11 (σ s 1j ) 2 + j=Nt j=1 1 2 j (w t lf ) j-1 w t 11 (σ t 1j ) 2 (10)

Model based bit allocation

As explained in [START_REF] Payan | An efficient bit allocation for compressing normal meshes with an error-driven quantization[END_REF], a bit allocation is generally processed with an iterative algorithm. At each iteration, we need to evaluate the distortion (and the associated bitrate) of each subband for one given set of quantization steps to check the convergence of the algorithm. In this way, one can achieve a real prequantization and then compute real distortion and bitrate for each subband. But such a method is time-consuming, particularly when dealing with large animations. One relevant method to overcome this problem is using statistical models to estimate distortion and bitrate of each subband without quantizing data. To perform such a model-based process, we must study the distribution of the processed wavelet coefficients. Figure 4 shows some results for the animated Venus.

We observe that the data follow a Gaussian distribution whatever the level of decomposition, and consequently we can use the model-based bit allocation proposed in [START_REF] Payan | An efficient bit allocation for compressing normal meshes with an error-driven quantization[END_REF].

Results

In this section we compare the proposed compression scheme to the following coders:

-The temporal wavelet coder proposed in [START_REF] Payan | Temporal wavelet-based geometry coder for 3d animations[END_REF] applied to the original irregular animations; -The previous spatio-temporal wavelet coder without an optimal bit allocation proposed in [START_REF] Payan | Remeshing and spatio-temporal wavelet filtering for 3d animations[END_REF].

To evaluate the quality of compressed animations, we use the RMSE metric based on Hausdorff distance given in [START_REF] Aspert | Mesh: Measuring errors between surfaces using the hausdorff distance[END_REF]. First, figures 5, 6, and 7 show the average reconstruction error (computed on all the frames) of these coders in function of different bitrates, for the animated Dolphin, Venus and Horse.

We observe that the proposed compression scheme based on a spatio-temporal decomposition and an optimal model-based bit allocation gives the best results. It clearly proves the interest of using a spatio-temporal decomposition and an optimal bit allocation process for such data instead of the previous allocation process proposed in [START_REF] Payan | Remeshing and spatio-temporal wavelet filtering for 3d animations[END_REF] In addition, figures 8, 9 and 10 show the evolution of the RMSE frame after frame when using the different coders, at given bitrates. We can observe that the proposed coder significantly improves the coding performances. With a smaller bitrate, the proposed coder gives a smaller reconstruction error for the majority of frames. Furthermore, it tends to produce a more regular quality over the frames (no peak of distortion, unlike the two other coders).

Conclusion and future works

In this paper, we have presented a new framework for compressing animations sharing the same connectivity at any frame. Inspired by works in geometry processing for static meshes, we have developed a spatio-temporal wavelet-based coding scheme for animated semi-regular meshes. Semi-regular animations offers the advantage of making the spatial wavelet-based analysis easier, since wavelets are more efficient on (semi-)regular grids than on irregular ones. The key contribution of this paper is the bit allocation process which optimizes the quantization of all the subbands of temporal and spatial wavelet coefficients created by the analysis tool. Furthermore, the statistics of the wavelet coefficients allows us to use a model-based algorithm, which makes the allocation process very fast. Experimentally, we show that the proposed method improves signifi-cantly the compression performances of the previous coders [START_REF] Payan | Temporal wavelet-based geometry coder for 3d animations[END_REF] and [START_REF] Payan | Remeshing and spatio-temporal wavelet filtering for 3d animations[END_REF], which makes our new coder particularly interesting. For future works we plan to improve the spatio-temporal filtering and add a motion estimation/compensation technique to further improve the performances of our coder.

Fig. 1 .

 1 Fig. 1. Overview of the complete proposed coding scheme.

  X-coordinates (second level of spatial decomposition).(e) Y-coordinates (second level of spatial decomposition).(f) Z-coordinates (second level of spatial decomposition).

Fig. 4 .Fig. 5 .

 45 Fig. 4. Typical distribution of the wavelet coefficients of the animated semi-regular Venus with the spatio-temporal decomposition used.

Fig. 6 .

 6 Fig. 6. Rate-Distortion curve for the animated Venus.

Fig. 7 .

 7 Fig. 7. Rate-Distortion curve for the animated Horse.

Fig. 8 .

 8 Fig. 8. RMSE-frame curve for the animated dolphin.

Fig. 9 .

 9 Fig. 9. RMSE-frame curve for the animated Venus.

Fig. 10 .

 10 Fig. 10. RMSE-frame curve for the animated Horse.

  One can obtain a multiresolution decomposition by subsequent filterings of the low frequency subbands. This decomposition finally consists in N t subbands of detail frames H t r t (where r t represents the temporal resolution), and the low frequency sequence L t Nt .

	h t 1,T /2 } and
	a low frequency subband L t 1 defined by T /2 approximation frames (temporal low frequency signal ): L t 1 = {l t 1,0 , l t 1,1 , ..., l t 1,T /2 }.

Acknowledgements

Original Venus Head and Horse models is courtesy of Cyberware. We are particularly grateful to J-H. Yang for providing us the original Dolphin model. The animated models Dolphin, Venus and Horse can be found at: http://www.i3s.unice.fr/˜fpayan/