Didac Barroso-Bergada

Charlie Pauvert

Jessica Vallance

Laurent Delière

David A Bohan

Marc Buée

Dr Corinne Vacher
email: corinne.vacher@inrae.fr

Microbial networks inferred from environmental DNA data for biomonitoring ecosystem change: strengths and pitfalls Running title (45 characters): Microbial networks for biomonitoring

Keywords: Environmental DNA, Metabarcoding, Community ecology, Ecosystem services, Microbial networks, Network inference, Network comparison

Environmental DNA contains information on the species interaction networks that support ecosystem functions and services. Next-Generation Biomonitoring proposes the use of this data to reconstruct ecological networks in real-time and then compute network-level properties to assess ecosystem change. We investigated the relevance of this proposal by assessing: (1) the replicability of DNA-based networks in the absence of ecosystem change; and, (2) the benefits and shortcomings of communityand network-level properties for monitoring change. We selected crop-associated microbial networks as a case study since they support disease regulation services in agroecosystems and analyzed their response to change in agricultural practice between organic and conventional systems. Using two statistical methods of network inference, we showed that network-level properties, especially βproperties, could detect change. Moreover, consensus networks revealed robust signals of interactions between the most abundant species, that differed between agricultural systems. These findings complemented those obtained with community-level data, that showed, in particular, a greater microbial diversity in the organic system. The limitations of network-level data included (i) the very high variability of network replicates within each system; (ii) the low number of network replicates per system, due to the large number of samples needed to build each network; and, (iii) the difficulty in interpreting links of inferred networks. Tools and frameworks developed over the last decade to infer and compare microbial networks are therefore relevant to biomonitoring, provided that the DNA metabarcoding datasets are large enough to build many network replicates and progress is made to increase network replicability and interpretation.

Introduction

Interactions among organisms and with their abiotic environment regulate the ecological processes underlying ecosystem services [START_REF] Mace | Biodiversity and ecosystem services: A multilayered relationship[END_REF]. Ecological interactions among organisms (e.g. predation, mutualism, parasitism) at a single point in space and time are usually represented as a network, with the organisms as nodes and the interactions as links [START_REF] Pocock | The robustness and restoration of a network of ecological networks[END_REF].

Current challenges focus on understanding how and why these networks vary in space and time [START_REF] Pellissier | Comparing species interaction networks along environmental gradients[END_REF][START_REF] Pilosof | The multilayer nature of ecological networks[END_REF], and which network properties should be conserved or enhanced to sustain ecosystem services [START_REF] Montoya | Emerging perspectives in the restoration of biodiversity-based ecosystem services[END_REF][START_REF] Raimundo | Adaptive Networks for Restoration Ecology[END_REF][START_REF] Tylianakis | Conservation of species interaction networks[END_REF]. Next-Generation Biomonitoring (NGB) proposes the reconstruction, automatically and in real-time, of ecological networks using the Next-Generation Sequencing (NGS) of environmental DNA (eDNA) data, and the analysis of network and community variation in space and time to detecting and explaining changes in ecosystem functions and services [START_REF] Baird | Biomonitoring 2.0: A new paradigm in ecosystem assessment made possible by next-generation DNA sequencing[END_REF][START_REF] Bohan | Next-Generation Global Biomonitoring: Large-scale, Automated Reconstruction of Ecological Networks[END_REF][START_REF] Derocles | Biomonitoring for the 21st Century: Integrating Next-Generation Sequencing Into Ecological Network Analysis[END_REF][START_REF] Makiola | Key Questions for Next-Generation Biomonitoring[END_REF]. However, before implementing NGB approaches on a large scale, we need more case studies demonstrating the utility of DNA-based networks and the meaning of their derived network metrics [START_REF] Compson | Network-Based Biomonitoring: Exploring Freshwater Food Webs With Stable Isotope Analysis and DNA Metabarcoding[END_REF]. The goal of the present study is to fill this gap.

NGB requires the reconstruction of replicated networks of ecological interactions as well as the development of statistical tools for their comparison and analysis. Theoretical frameworks have been developed for the comparison of ecological networks between contrasted environmental conditions or along environmental gradients [START_REF] Delmas | Analysing ecological networks of species interactions[END_REF][START_REF] Pellissier | Comparing species interaction networks along environmental gradients[END_REF][START_REF] Poisot | The dissimilarity of species interaction networks[END_REF][START_REF] Tylianakis | Ecological Networks Across Environmental Gradients[END_REF]. By analogy with the α-and β-diversity of ecological communities, these frameworks define α-and β-properties for ecological networks as whole-network metrics (e.g. connectance) and dissimilarities between pairs of networks, respectively [START_REF] Pellissier | Comparing species interaction networks along environmental gradients[END_REF]. Community-and network-level metrics can be used to assess the impact of environmental changes on the number, identity and abundance of the species forming ecological communities, and on the structure, type and strength of their interactions, respectively. They have for instance been used to evaluate the impact of agricultural practices [START_REF] Morriën | Soil networks become more connected and take up more carbon as nature restoration progresses[END_REF] , that are a key driver of global change [START_REF] Tilman | Agricultural sustainability and intensive production practices[END_REF], on species diversity [START_REF] Tuck | Land-use intensity and the effects of organic farming on biodiversity: A hierarchical metaanalysis[END_REF] and on pest and disease regulation services supported by species interactions [START_REF] Ma | Ecological networks reveal resilience of agro-ecosystems to changes in farming management[END_REF][START_REF] Macfadyen | Parasitoid control of aphids in organic and conventional farming systems[END_REF][START_REF] Tylianakis | Habitat modification alters the structure of tropical host-parasitoid food webs[END_REF].

Networks of interactions among microorganisms appear as suitable tools for NGB for at least three reasons: NGS techniques are the current rule for studying microbial communities [START_REF] Bálint | Millions of reads, thousands of taxa: Microbial community structure and associations analyzed via marker genesa[END_REF]; microorganisms are present in all Earth ecosystems; and, microbial interactions are crucial to ecosystem functioning, human life and well-being [START_REF] Gilbert | Life in a World without Microbes[END_REF][START_REF] Zhu | Changes in the environmental microbiome in the Anthropocene[END_REF].

Network ecology, which originates from the study of trophic links between macroorganisms [START_REF] Ings | Ecological networks -Beyond food webs[END_REF], initially ignored interactions with and among smaller organisms [START_REF] Lafferty | Parasites dominate food web links[END_REF]. But, with increasing evidence of the contribution of microbial interactions to biogeochemical cycles [START_REF] Falkowski | The microbial engines that drive earth's biogeochemical cycles[END_REF], plant diversity and productivity (van der Heijden, [START_REF] Van Der Heijden | The unseen majority: Soil microbes as drivers of plant diversity and productivity in terrestrial ecosystems[END_REF], and disease regulation in soils [START_REF] Berendsen | The rhizosphere microbiome and plant health[END_REF], plants and animals [START_REF] Brader | Ecology and Genomic Insights into Plant-Pathogenic and Plant-Nonpathogenic Endophytes[END_REF][START_REF] Hacquard | Interplay Between Innate Immunity and the Plant Microbiota[END_REF][START_REF] Vayssier-Taussat | Shifting the paradigm from pathogens to pathobiome: new concepts in the light of meta-omics[END_REF], microbial networks are now considered key to the understanding of ecosystem functioning [START_REF] De Vries | Soil bacterial networks are less stable under drought than fungal networks[END_REF][START_REF] Karimi | Microbial diversity and ecological networks as indicators of environmental quality[END_REF][START_REF] Wagg | Fungal-bacterial diversity and microbiome complexity predict ecosystem functioning[END_REF]. However, given that microbial networks inferred from eDNA data only represent hypothesized interactions among microbial species but rather statistical associations among molecular units that only represent putative signals for microbial interactions [START_REF] Faust | Microbial interactions: From networks to models[END_REF][START_REF] Röttjers | From hairballs to hypotheses-biological insights from microbial networks[END_REF][START_REF] Vacher | The Phyllosphere: Microbial Jungle at the Plant-Climate Interface[END_REF], it is crucial to evaluate the relevance of the derived network properties to the assessment of change in ecosystem functioning.

In this study, we analyzed the relevance of microbial network properties to NGB, by assessing (1) the replicability of microbial networks inferred from eDNA data in the absence of ecosystem change, and (2) the benefits and shortcomings of community-level and network-level properties for detecting change. We focused on crop-associated microbial networks since they support disease regulation services in agroecosystems [START_REF] Toju | Core microbiomes for sustainable agroecosystems[END_REF], and analyzed their response to change in agricultural practice (conventional vs organic farming). We inferred microbial networks from eDNA sampled from replicated agricultural plots by using two classical methods of network inference, SparCC [START_REF] Friedman | Inferring Correlation Networks from Genomic Survey Data[END_REF] and SPIEC-EASI [START_REF] Kurtz | Sparse and Compositionally Robust Inference of Microbial Ecological Networks[END_REF]. We then computed α-and β-diversity metrics at the community-and network-level to identify the level that best captures change in agricultural practice, by using grapevine and its foliar microorganisms as the case study. These results are then used to discuss those tools and frameworks that are best adapted to NGB approaches.

Materials and methods

Study site and sampling design

Grapevine leaf samples were collected on September 10, 2015, from an experimental vineyard (Figure 1) located near Bordeaux (INRA, Villenave d'Ornon, France; 44°47'32.2"N 0°34'36.9"W). The experimental vineyard was planted in 2011 and was designed to compare three cropping systems: sustainable conventional agriculture (CONV), organic farming (ORGA) and pesticide-free farming (RESI) [START_REF] Delière | Le réseau Ecoviti Bordeaux expérimente des systèmes de culture viticoles « bas-intrants[END_REF]. The Vitis vinifera L. cultivar Merlot noir grafted onto a 3309 C rootstock was used in both the CONV and ORGA cropping systems. Only the CONV and ORGA systems, that used the same cultivar but different phytosanitary treatments, were compared in the present study to avoid multiplying the sources of variation between systems. RESI used a resistant cultivar, which has 2 Quantitative Trait Loci of partial resistance to downy mildew and total resistance to powdery mildew.

The experiment had a randomized block design [START_REF] Schielzeth | Nested by design: Model fitting and interpretation in a mixed model era[END_REF] consisting of three blocks, each composed of three plots, one for each of the cropping systems tested. Each plot covered an area of 2100 m² and was composed of 20 rows of 68 vines each, with 1.60 m between rows and 0.95 m between vines in a single row.

CONV plots were managed according to the general principles of integrated pest management (IPM), as listed in Appendix III of the 2009/128/EC Directive (European Commission, 2009). ORGA plots were managed according to European Council Regulation (EC) No 834/2007(European Council, 2007).

ORGA plots were treated with copper and sulfur-based products, whereas additional phytosanitary products were allowed in CONV plots (Table S1). The cropping systems differed in terms of the types of pesticides applied and the timing of applications, but not in terms of doses (Table S1). All products and active ingredients were applied between the end of April and mid-August of 2015. Grapes were harvested on September 10, 2015. The disease incidence and severity at harvest were higher in CONV plots than in ORGA plots for both powdery mildew (caused by the fungal pathogen Erysiphe necator) and black rot (caused by the fungal pathogen Guignardia bidwellii). Downy mildew symptoms (caused by the oomycete pathogen Plasmopara viticola) did not differ significantly between the cropping systems (Table S2).

Grapevine leaves were collected in the two hours prior to grape harvest, from 20 vines per plot in the CONV and ORGA plots (Figure 1). We attempted to avoid edge effects by selecting the 20 vines from the center of each plot. The third leaf above the grapes was collected from each vine, placed in an individual bag and immediately transported to the laboratory. In total, 120 leaves, corresponding to 1 leaf × 20 vines × 3 plots × 2 cropping systems, were collected. Leaves were processed on the day of collection, with sterilized tools in the sterile field of a MICROBIO electric burner (MSEI, France).

Three contiguous discs of 6 mm diameter were cut from the center of each leaf, approximately 2 cm from the midrib. They were placed in the well of a sterile DNA extraction plate. The leaf disks were then freeze-dried overnight (Alpha 1-4 DA Plus, Bioblock Scientific).

DNA extraction and sequencing

Leaf disks (Figure 1) were ground with a single-glass ball mill (TissueLyser II, Qiagen) and DNA was then extracted with a CTAB chloroform/isoamyl alcohol (24:1) protocol. A dozen "empty" wells (i.e. containing nothing but extraction reagents) were included on each plate as negative control samples for DNA extraction. Three of these negative control samples were randomly selected and pooled before sequencing. Three replicates of a fungal mock community, each consisting of an equimolar pool of DNA from 189 pure fungal strains, were also included as positive control samples [START_REF] Pauvert | Bioinformatics matters: The accuracy of plant and soil fungal community data is highly dependent on the metabarcoding pipeline[END_REF].

The nuclear ribosomal internal transcribed spacer (ITS) region, which is considered to be the universal barcode region for fungi [START_REF] Schoch | Nuclear ribosomal internal transcribed spacer (ITS) region as a universal DNA barcode marker for Fungi[END_REF], was then amplified with the ITS1F (5'-CTTGGTCATTTAGAGGAAGTAA-3', [START_REF] Gardes | ITS primers with enhanced specificity for basidiomycetes -application to the identification of mycorrhizae and rusts[END_REF] and ITS2 (5'-GCTGCGTTCTTCATCGATGC-3', [START_REF] White | Amplification and Direct Sequencing of Fungal Ribosomal Rna Genes for Phylogenetics[END_REF] primer pair, which targets the ITS1 region. PCR was performed in an Eppendorf thermocycler (Eppendorf), with a reaction mixture (25 µl final volume) consisting of 0.04 U Taq polymerase (SilverStar DNA polymerase, Eurogentec), 1X buffer, 2 mM MgCl2, 200 µM of each dNTP, 0.2 µM of each primer, 1 ng.µl -1 bovine serum albumin (New England BioLabs) and 2 µl DNA template. A pseudo-nested PCR protocol was used, with the following cycling parameters: enzyme activation at 95°C for 2 min; 20 (1 st PCR with regular primers; Table S3) and then 15 (2 nd nested PCR with pre-tagged primers; Table S3) cycles of denaturation at 95°C for 30 s, 53°C for 30 s, 72°C for 45 s; and a final extension phase at 72°C for 10 min. "Empty" wells (i.e. containing nothing but PCR reagents) were included on each plate as a negative control for PCR. Three negative control samples were randomly selected and pooled before sequencing. In addition, the PCR product of one sample per plot was split in two, with each half of the sample sequenced independently to serve as technical replicates for sequencing, hence forming six pairs of technical replicates (one per plot).

We checked the quality of all the PCR products by electrophoresis in 2% agarose gels. A total of 123 samples were sent to sequencing, corresponding to 112 well-amplified leaf samples, 6 technical replicates, 1 pooled negative extraction control, 1 pooled negative PCR control and 3 mock community replicates. PCR products were purified (CleanPCR, MokaScience), multiplex identifiers and sequencing adapters were added, and library sequencing on an Illumina MiSeq platform (v3 chemistry, 2×250 bp) and sequence demultiplexing (with exact index search) were performed at the Get-PlaGe sequencing facility (Toulouse, France).

Bioinformatic analysis

Based on the mock community included in the sequencing run, we found that analyzing single forward (R1) sequences with DADA2 [START_REF] Callahan | DADA2: High-resolution sample inference from Illumina amplicon data[END_REF] was a good option for fungal community characterization [START_REF] Pauvert | Bioinformatics matters: The accuracy of plant and soil fungal community data is highly dependent on the metabarcoding pipeline[END_REF]. This pipeline fully exploits the resolution of molecular barcodes [START_REF] Callahan | DADA2: High-resolution sample inference from Illumina amplicon data[END_REF], which is a desired feature in microbial network inference. Indeed, the taxonomic resolution of the nodes should be fine enough to discern the variation in ecological interactions between microbial strains [START_REF] Röttjers | From hairballs to hypotheses-biological insights from microbial networks[END_REF]. Using DADA2 v1.6, we retained only R1 reads with less than one expected error (based on quality scores; [START_REF] Edgar | Error filtering, pair assembly and error correction for next-generation sequencing reads[END_REF] that were longer than 100 bp, and we then inferred amplicon sequence variants (ASV) for each sample. Chimeric sequences were identified by the consensus method of the removeBimeras function. Taxonomic assignments were performed with RDP classifier [START_REF] Wang | Naïve Bayesian classifier for rapid assignment of rRNA sequences into the new bacterial taxonomy[END_REF], implemented in DADA2 and trained with the UNITE database v. 7.2 (UNITE Community 2017). Only ASVs assigned to a fungal phylum were retained. The ASV table was then filtered as described by [START_REF] Galan | 16S rRNA Amplicon Sequencing for Epidemiological Surveys of Bacteria in Wildlife[END_REF] with a custom script (https://gist.github.com/cpauvert/1ba6a97b01ea6cde4398a8d531fa62f9) that removed ASVs from all samples for which the number of sequences was below the crosscontamination threshold, defined as their maximum number in negative control samples. Finally, we checked the compositional similarity of the six pairs of technical replicates, in terms of both ASV occurrence and relative abundance (Figure S1), and we removed for each pair of technical replicates the replicate with the lowest number of sequences. We also removed the controls. Therefore, the final ASV table contained 1116 ASVs, 112 leaf samples and 4,760,068 high-quality sequences.

Statistical analyses

Statistical analyses were performed with R software v3.4.1 (R Core Team, 2018), with the packages lme4 v1.1-19 [START_REF] Bates | Fitting linear mixed-effects models using lme4[END_REF], vegan v2.5-5 (Oksanen et al., 2019), permute v0.9-5 [START_REF] Simpson | permute: Functions for Generating Restricted Permutations of Data (0.9-5). R package[END_REF]), phyloseq v1.24.2 (McMurdie & Holmes, 2013) including the DESeq2 extension v1.20.0 [START_REF] Love | Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[END_REF]), NST v2.0.4 (Ning et al., 2019), and igraph v1.2.4.1 [START_REF] Csardi | The igraph software package for complex network research[END_REF]. Data were manipulated and plots were created with reshape2 v1.4.3, plyr v1.8.4 and ggplot2 v3.2.0 [START_REF] Wickham | ggplot2 -Elegant Graphics for Data Analysis[END_REF], cowplot v0.9.4 [START_REF] Wilke | cowplot: Streamlined Plot Theme and Plot Annotations for "ggplot2[END_REF], ggraph v1.0.2 [START_REF] Pedersen | ggraph: An Implementation of Grammar of Graphics for Graphs and Networks[END_REF] and VennDiagram v1.6.20 [START_REF] Chen | VennDiagram: Generate High-Resolution Venn and Euler Plots[END_REF].

Effect of cropping system on community α-diversity

Three community α-diversity properties were computed for each sample: richness, diversity and evenness of fungal communities (Table 1). Generalized linear mixed models (GLMMs) were then used to test the effect of the cropping system on these properties. The models included the cropping system as a fixed treatment effect and the sampling depth (defined as the total number of raw sequences per sample) as an offset [START_REF] Bálint | Relocation, high-latitude warming and host genetic identity shape the foliar fungal microbiome of poplars[END_REF][START_REF] Mcmurdie | Waste Not, Want Not: Why Rarefying Microbiome Data Is Inadmissible[END_REF]. For every property, we compared the likelihood of a full model, including the block and its interaction with the cropping system as random effects and a simplified model, including only the block factor as a random effect. Community richness was defined as the number of ASVs per sample. We used a logarithmic link function to model these count data, assuming a negative binomial distribution to deal with overdispersion [START_REF] Zuur | Mixed Effects Modelling for Nested Data[END_REF]. Community diversity was measured with the Inverse Simpson index [START_REF] Simpson | Measurement of diversity[END_REF] and modeled with a Gaussian distribution and the logarithmic link function. Evenness was estimated with Pielou's index [START_REF] Pielou | The measurement of diversity in different types of biological collections[END_REF] and modeled with a Gaussian distribution and the logarithmic link function. The offset was transformed according to the link function. The significance of the fixed treatment effect was finally assessed with the Wald χ² test [START_REF] Bolker | Generalized linear mixed models: a practical guide for ecology and evolution[END_REF]. Moreover, to investigate whether foliar fungal pathogens of grapevine were responsible for variations in community α-diversity properties, we fitted the models by including the relative abundance of sequences assigned to the Erysiphe genus (which includes Erysiphe necator, the causal agent of powdery mildew; [START_REF] Armijo | Grapevine pathogenic microorganisms: Understanding infection strategies and host response scenarios[END_REF] and the Guignardia genus (which includes Guignardia bidwelli, the causal agent of black rot) as fixed additive effects.

Effect of cropping system on community β-diversity

Two community β-diversity properties were calculated for each pair of samples: the quantitative Jaccard dissimilarity and the binary Jaccard dissimilarity (Table 1). Permutational analyses of variance (PERMANOVAs; [START_REF] Anderson | Permutation tests for univariate or multivariate analysis of variance and regression[END_REF] were then used to evaluate the effect of the cropping system on these compositional dissimilarities. The models included cropping system, sampling depth (logtransformed), block and their interaction as fixed effects. ASVs differing in abundance between cropping systems were identified with DESeq2 [START_REF] Love | Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[END_REF], by calculating the likelihood ratio between a full model including block and cropping system as fixed effects and a simplified model including only the block factor. The estimated fold-changes in abundance were considered significant if the p-value was below 0.05 after Benjamini and Hochberg adjustment. Moreover, to understand better the processes shaping community structure, the relative contribution of deterministic and stochastic processes in community assembly was assessed by following the framework defined by [START_REF] Ning | A general framework for quantitatively assessing ecological stochasticity[END_REF]. This method provides statistics for each sample, named the Normalized Stochasticity Ratio (NST), that ranges from 0 to 100, where 0 means a completely deterministic assembly process and 100 a completely stochastic assembly process. NST was calculated using the tNST function with the binary and quantitative Jaccard dissimilarity indices, the FE null model, and other parameters by default values. We used the FE null model (SIM2 in [START_REF] Gotelli | Null model analysis of species co-occurrence patterns[END_REF]) because it is the most appropriate for comparing standardized samples that have been collected in areas of homogeneous habitat, such as vineyards. This null model reshuffles ASV occurrences among samples by considering that all samples are equally probable. NST values were calculated for each cropping system and then compared using permutational analysis of variance with the nst.panova function.

Network inference

Fungal association networks were inferred at plot level (Figure 1) with two widely-used methods of microbial network inference: the SparCC algorithm [START_REF] Friedman | Inferring Correlation Networks from Genomic Survey Data[END_REF] implemented in FastSpar [START_REF] Watts | FastSpar: Rapid and scalable correlation estimation for compositional data[END_REF] with default SparCC values; and, the SPIEC-EASI method [START_REF] Kurtz | Sparse and Compositionally Robust Inference of Microbial Ecological Networks[END_REF] using the MB procedure of edge selection. Both methods try to deal with the compositional nature of metabarcoding data. In a metabarcoding dataset, the total number of sequences per sample is arbitrary, imposed by the sequencer. Sequence counts contain only relative abundance information for species. Methods that do not take this feature into account can result in the identification of artifactual associations [START_REF] Gloor | Microbiome datasets are compositional: And this is not optional[END_REF].

Both SparCC [START_REF] Friedman | Inferring Correlation Networks from Genomic Survey Data[END_REF] and SPIEC-EASI [START_REF] Kurtz | Sparse and Compositionally Robust Inference of Microbial Ecological Networks[END_REF] attempt to overcome this bias using log ratios of counts. For each method of network inference, ten networks per plot were constructed by varying the percentage P of the ASVs included in the network (with P ranging from 10% to 100% of the most abundant ASVs in the plot). We varied P because we expected that it would influence the replicability of the networks. We expected, in particular, the networks built from only the most abundant ASVs to be more replicable. For the same reason, networks were also inferred after aggregating ASVs at the genus level and removing ASVs that were not taxonomically assigned at this level. In all cases, the inferred microbial networks had ASVs as nodes and a positive or negative link between ASVs in cases of significant associations between abundance.

Effect of cropping system on network α-properties

Six network α-properties were calculated for each inferred network: number of links, network density, number of connected components, diameter of the largest component, mean node degree and proportion of negative links (Table 1). The effect of the cropping system on these properties was investigated by performing Wilcoxon rank-sum tests for all values of P. The Benjamini-Hochberg procedure was used to correct p-values for multiple testing.

Effect of cropping system on network β-properties

Four network β-properties were calculated for each pair of inferred networks (Table 1). The topological distance between networks was calculated with the D index defined by [START_REF] Schieber | Quantification of network structural dissimilarities[END_REF].

Schieber's D, when applied to binary networks (i.e. with unweighted links) captures global and local structural dissimilarities between networks, by comparing node connectivity patterns across scales.

The dissimilarity of associations between networks, βWN, according to the framework described by Poisot et al. (2012), was then calculated for all pairs of networks with the binary Jaccard dissimilarity index. βWN was then partitioned into two components (Poisot et al., 2012): the dissimilarity of associations between ASVs common to both networks (βOS) and the dissimilarity of associations due to species turnover (βST). In contrast to the Shieber's D index that evaluates how nodes are connected to neighboring nodes and to more distant nodes, these three metrics compare lists of pairwise associations between nodes. PERMANOVA was used to evaluate the effect of the cropping system on the topological distance between networks (D) and the dissimilarity of associations between networks (βWN, βST and βOS). The models included cropping system, the percentage of ASV, P, and their interactions as fixed effects. The permutations (n=999) were constrained within blocks. Finally, for each network inference method and every value of P, consensus networks containing only the shared associations between the three network replicates within a cropping system were built to identify robust associations that could indicate ecological interactions between fungal strains. The number of shared associations between the three network replicates were compared to those obtained between three random networks simulated with the same nodes and the same number of links. The significance of shared associations was evaluated with a pseudo p-value, estimated from 999 simulations and defined as the probability that the three random networks shared more associations than the three inferred networks [START_REF] Morlon | Effects of trophic similarity on community composition[END_REF].

Results

Among the 15 community-and network-level properties computed (Table 1), 7 indicated differences between the organic (ORGA) and the conventional (CONV) system.

All community α-properties detected system change

All three community α-diversity properties -richness, diversity and evenness (Table 1) -were significantly higher in ORGA than CONV plots (Figure 2 A-C and Table S4).Community richness, for example, equaled on average 39.69 fungal ASVs per sample in ORGA plots vs 36.40 in CONV plots, with each sample representing 0.85 cm2 of a single leaf tissue. Including the interaction between the cropping system and the block did not significantly increase the likelihood reported by the GLMM models, indicating that changes in community α-diversity properties due to the cropping system were consistent across blocks.

In contrast to our expectations, none of the community α-diversity properties was influenced by pathogen relative abundance (Table S5). Pathogen abundance within each sample was estimated as the proportion of sequences assigned to the Erysiphe genus and ranged between 0% and 36.34%, with an average of 1.12% per sample. No ASV was assigned to the Guignardia genus and this variable was therefore not included in the models.

All community β-properties detected system change

The two community β-diversity properties analyzed in this study -the quantitative and binary Jaccard indices (Table 1) -detected significant differences in community composition between systems (Table 2). The cropping system was a major driver of both ASV relative abundance (Figure 2D) and ASV presence-absence (Figure 2E), as indicated by the quantitative and binary Jaccard indices, respectively. It explained 7.6% of the variance in ASV relative abundance and 4.5% of the variance in ASV presenceabsence (Table 2). The block effect was also significant, indicating that there were spatial variations in community composition at the scale of the experiment. The block explained 4.3% of the variance in ASV relative abundance, and 2.6% of the variance in ASV presence-absence (Table 2). There were also large differences in composition among samples within a plot, as indicated by the high percentage of unexplained variance (78.2% for the quantitative Jaccard index and 85.7% for the binary Jaccard index) (Table 2).

In line with these results, we found that the stochasticity in ASV presence-absence was very high in both the ORGA and CONV systems (NST=78.4% and 94.8%, respectively). Nevertheless, it decreased markedly when the relative abundance of ASVs (NST=29.3% and 33.6%, respectively) was taken into account (Table S6), probably because the ASV, assigned to Aureobasidium sp. (Table 3) was the most abundant, represented more than half of the total number of sequences and was highly abundant in all samples. Stochasticity in ASV presence-absence was significantly higher in CONV plots (Table S6). A similar trend, although non significant, was observed for ASV relative abundance, suggesting that communities in ORGA plots were more stable, in addition to being richer (Figure 2A).

Overall, the foliar fungal communities were dominated by Ascomycota in both ORGA (87.2% of sequences) and CONV (96.8%) plots. About one-fourth of ASVs (249 over 1116) were shared between cropping systems. These shared ASVs were the most abundant, representing 98.97% of the total number of sequences. Fourteen ASVs differed significantly in abundance between the cropping systems according to differential abundance analysis performed with DeSeq2 (Figure 2F). For instance, the causal agent of grapevine powdery mildew, Erysiphe necator, which was among the 10 most abundant fungal species, was significantly more abundant in CONV than in ORGA plots (Figure 2F), according to both the visual records of disease symptoms (Table S2) and metabarcoding data (2% versus less than 0.1%; Table 3). The highest abundance of this major grapevine pathogen in samples of CONV plots was however not responsible for their lower α-diversity (Figure 2 A-C and Table S5).

Differential abundance analysis also revealed that three other ASVs were significantly more abundant in CONV plots, whereas 10 other ASVs, including several yeast species (from the genera Vishniacozyma, Sporobolomyces and Filobasidium), were significantly more abundant in ORGA plots (Figure 2F).

None of the network α-properties detected system change

For each method of network inference, we obtained sixty fungal association networks (SparCC: Figure 3A; SPIEC-EASI: Figure S2A), each corresponding to one of the six vineyard plots (Figure 1) and one of the ten values of the percentage P of most abundant ASVs included in the network. Whatever the network inference method, none of the six network α-properties (Table 1) differed between cropping systems (Table S7 andS8), but all were significantly correlated with P (Table S9 andS10).

Four network α-properties had consistent variations with P between the two methods: the total number of links (L), the number of connected components (CC), the network connectance (C) and the average node degree (DEG) (Tables S9 andS10). Increasing the number of ASVs included in the network increased the total number of links, linked the connected components (hence reducing CC) and increased the average node degree. This consistent increase in average node degree with P, however, masked some differences between methods. With SPIEC-EASI, node degree increased more in abundant ASVs, yielding a significant, positive relationship between ASV relative abundance and node degree at P=100% (Figure S3). This was not the case in SparCC (Figure S3). Despite this difference, the network connectance decreased with both methods of network inference, consistent with their sparsity assumption [START_REF] Friedman | Inferring Correlation Networks from Genomic Survey Data[END_REF][START_REF] Kurtz | Sparse and Compositionally Robust Inference of Microbial Ecological Networks[END_REF].

Half of the network β-properties detect system change

Only two network β-properties, of the four computed (Table 1), differed significantly between cropping systems whatever the network inference method. As with the network α-properties , the topological dissimilarity between networks, measured with the Shieber's D index [START_REF] Schieber | Quantification of network structural dissimilarities[END_REF], did not differ between cropping systems but was influenced by P, irrespective of the network inference method (Table 4 andTable S11). These results are consistent with the results obtained for node degree and network connectance, which are components of the D index and also vary with P but do not differ between cropping systems (Tables S7 to S10).

By contrast, cropping system had a significant effect on the overall dissimilarity of associations (βWN) and the dissimilarity of associations between shared ASVs (βOS) for both SparCC networks (Table 4 and Figure 3C) and SPIEC-EASI networks (Table S11 and Figure S2C). Cropping system also had a significant effect on the dissimilarity of associations due to ASV turnover (βST), but only in SparCC networks and only in interaction with P (Table 4 andTable S11). These findings suggest that network variation between cropping systems are due to the turnover in associations (captured by βOS), rather than the turnover in ASVs (captured by βST), and show that the network β-properties defined in the theoretical ecology framework by Poisot et al. (2012) can be used to detect differences between cropping systems.

Network replicates within each system were highly variable but shared links

Network replicates varied considerably within a cropping system, whatever the network inference method (Figures 4 andS4). When all ASVs were used for network construction with SparCC (P=100%), only 3 associations were common to all three network replicates of the ORGA system, although 80 ASVs were shared between the three network replicates (Figure 4). Only 5 were common to all three network replicates of the CONV system, although 81 ASVs were shared between the three network replicates (Figure 4). Similar results were obtained with SPIEC-EASI, with 1 and 5 shared associations, respectively (Figure S4).

High variability of network replicates within a cropping system was observed for all values of P and was not reduced by the aggregation of ASVs at the genus level nor by the consideration of only the most abundant ASVs, in contrast with our expectation. The networks inferred from only the most abundant ASV or the most abundant genera (P=10% or 20%) barely shared any associations (Table 5). These low numbers of shared associations between the three networks replicates were, however, generally significantly higher than expected from three random networks having the same number of nodes and links (Table 5). The number of shared associations between the three network replicates ranged between 0 and 7, depending on P, and the network inference method (Table 5), while the average number of shared associations between the random networks ranged between 0 and 1.2, suggesting that consensus networks within a cropping system (Figure S5) do contain robust associations but these are few in number.

Five of nine consensus associations were also found by both methods of network inference. The SparcCC and SPIEC-EASI consensus networks obtained for P=100% in the ORGA system shared a negative association between the dominant ASV, assigned to Aureobasidium sp., and the third most abundant ASV in the ORGA system, assigned to Cladosporium ramotenellum (Figure S5 and Table 3), as an example The consensus networks obtained for the CONV system also shared a negative association between the dominant ASV, assigned to Aureobasidium sp., and the third most abundant ASV in the CONV system, assigned to Epicoccum nigrum (Figure S5 and Table 3). Three positive associations were also shared by the SparCC and SPIEC-EASI consensus networks in the CONV system (Figure S5). No association was shared between the two cropping systems, whatever the network inference method and despite 44 ASVs being shared (Figures 4 andS4 ), confirming the significant turnover in associations detected with βOS (Table 4 andTable S11).

Discussion

The functioning of ecosystems, like that of all complex systems, emerges from the interaction links between its components, and cannot be deduced from a simple listing of organisms [START_REF] Newman | The Structure and Dynamics of Networks[END_REF]. The concept of Next-Generation Biomonitoring (NGB) builds on this property of complex systems and proposes the use of networks of species interactions, rather than a simple list of species, to monitor changes in ecosystem functioning. It also proposes that this could be done via the automatic reconstruction of ecological networks from DNA metabarcoding data [START_REF] Bohan | Next-Generation Global Biomonitoring: Large-scale, Automated Reconstruction of Ecological Networks[END_REF]. In the present study, we focused on microbial association networks as a tool for ecosystem monitoring because microbial networks are present in all ecosystems, contribute to ecosystem functioning, and many methods exist to reconstruct them from DNA metabarcoding data [START_REF] Weiss | Correlation detection strategies in microbial data sets vary widely in sensitivity and precision[END_REF][START_REF] Dohlman | Mapping the microbial interactome: Statistical and experimental approaches for microbiome network inference[END_REF]. We assessed the relevance of microbial networks for NGB approaches using two criteria: (1) their replicability in the absence of environmental change; and,

(2) their ability to better detect environmental change than properties at the microbial community level. We focused on a major driver of environmental change, agricultural practices (conventional versus organic agriculture). Our results demonstrated that: (1) microbial network replicates were highly variable within the same set of environmental conditions; and, ( 2) some network-level metrics, but not all, could detect environmental change. By contrast, all community-level metrics revealed clear-cut changes in the microbial communities in response to environmental change (Table 1).

The high variability of network replicates within an environmental condition (i.e., in our study, the same cropping system) is the most surprising result of our study. When the whole metabarcoding dataset was used to build the networks, each network replicate was composed of about 160 nodes (fungal ASVs, in our study) and about 3500 links between these nodes (corresponding to co-occurrence or co-exclusion relationships between these fungi). The three network replicates shared half of their nodes but less than 5 links (Figure 4). Four non-mutually exclusive hypotheses can be put forward to explain this result. First, the variability in microbial associations may reflect real ecological variability.

Different assemblages of fungal taxa could play the same role in the ecosystem because of the functional redundancy of the taxa [START_REF] Louca | Decoupling function and taxonomy in the global ocean microbiome[END_REF]. There would thus be several assemblages, involving different associations of taxa (and thus different networks), adapted to the same cropping system. Second, the relative abundances of fungal taxa, from which the networks are built, could vary within the same environmental condition because of ecological stochasticity. The fungal communities were, like most ecological communities [START_REF] Mcintosh | Raunkiaer's "Law of Frequency[END_REF], composed of a small number of ubiquitous species and a large number of rare species whose presence varied greatly, probably because of the large degree of stochasticity in the deposition of fungal spores (Peay & Bruns, 2014). This high stochasticity in the composition of the rare microbiome may be responsible for the large number of associations that are unique to each network replicate and explain why the few shared associations involved abundant taxa. Third, the relative abundances of fungal taxa, upon which the networks are built, could vary within the same environmental condition because of methodological biases.

Distortions in taxa abundance may be generated at each step of the DNA metabarcoding process, from the collection of samples to their sequencing, and at each step of the bioinformatic processing of the sequences [START_REF] Ruppert | Past, present, and future perspectives of environmental DNA (eDNA) metabarcoding: A systematic review in methods, monitoring, and applications of global eDNA[END_REF]. The fungal ITS region, which was used here as a barcode [START_REF] Schoch | Nuclear ribosomal internal transcribed spacer (ITS) region as a universal DNA barcode marker for Fungi[END_REF], is highly variable in terms of length, sequence and number of copies [START_REF] Nilsson | Intraspecific ITS variability in the Kingdom Fungi as expressed in the international sequence databases and its implications for molecular species identification[END_REF][START_REF] Lofgren | Genome-based estimates of fungal rDNA copy number variation across phylogenetic scales and ecological lifestyles[END_REF], and these features could have increased the variability in the sequence data. Metabarcoding data are inherently noisy and this noise may explain why many associations are unique to a network replicate. Fourth, environmental conditions, which we consider homogeneous within a culture system, may not be homogeneous for microorganisms. Our experimental system and sampling protocols were designed to limit environmental variations within a cropping system. The vineyard plots were adjacent to each other and planted with grapevine clones. Moreover, we collected all leaves in less than two hours and controlled for the position of the sampled leaf on the vine.

Nonetheless, the significant block effects in community composition indicate that the fungal 20 communities were spatially structured at the scale of the experiment, which could account for spatial variations in networks. This poses a fundamental problem for biomonitoring approaches. The changes we want to monitor, which are generally large-scale changes in ecosystem functioning induced by human activities, may not necessarily be those to which microbial communities and networks respond.

Our study also highlighted a major pitfall of network comparison analyses, which is the lack of statistical power due to the low number of network replicates. To evaluate the effect of the cropping system, we had 56 replicates per system at the community level, but only 3 at the network level.

Indeed, several communities are needed to build a single network. This could explain why all community-level α-properties, but no network-level α-property, detected changes triggered by the cropping system. Despite this lack of statistical power, β-properties of microbial networks differed significantly between cropping systems, revealing a difference in microbial associations between organic and conventional systems. These differences were significant when network pairwise comparisons were based on shared taxa only, suggesting that the differences between organic and conventional networks were not only due to the turnover of taxa between cropping systems, but to re-associations of taxa. Overall, these results show that microbial networks inferred from DNA metabarcoding data can be used to detect changes in ecosystems if they are analyzed with network comparison tools defined by theoretical ecology (Pelissier et al., 2018;Poisot et al. 2012). They also suggest that β-properties of networks are better indicators of change than α-properties.

Our study also allowed us to compare two microbial network inference methods, SparCC [START_REF] Friedman | Inferring Correlation Networks from Genomic Survey Data[END_REF] and SPIEC-EASI [START_REF] Kurtz | Sparse and Compositionally Robust Inference of Microbial Ecological Networks[END_REF]. The results obtained with the two methods were, overall, encouragingly consistent. The variability of network replicates within a culture system was very high, regardless of the inference method used. The number of associations per network was lower with SPIEC-EASI than with SparCC (about 800 vs. 3500), probably because SPIEC-EASI infers partial correlations, discarding the indirect associations retained by SparCC [START_REF] Kurtz | Sparse and Compositionally Robust Inference of Microbial Ecological Networks[END_REF]. However, the number of associations shared between network replicates was very low in both cases (less than 5) (Figure 4 and S4). SPIEC-EASI found slightly fewer shared associations than SparCC, especially when the number of nodes were reduced by filtering on taxa abundance or taxonomic aggregation.

However, the shared associations detected by SPIEC-EASI had a higher level of significance (Table 5).

These results are in line with previous benchmarking studies showing the lower performance of SparCC compared to other methods of network inference, including SPIEC-EASI (Rötjjers & Faust, 2018;[START_REF] Hirano | Difficulty in inferring microbial community structure based on co-occurrence network approaches[END_REF] even though SparCC seems to work in low diversity communities [START_REF] Weiss | Correlation detection strategies in microbial data sets vary widely in sensitivity and precision[END_REF]. Both methods, however, revealed very similar consensus associations within each cropping system. Nine associations, in total, were shared by the network replicates and 5 were found by both methods. Although they involved ubiquitous fungal species that have been frequently detected on grapevine, such as Aureobasidium pullulans, Epicoccum nigrum and Cladosporium ramotenellum [START_REF] Martini | DNA-dependent detection of the grapevine fungal endophytes aureobasidium pullulans and epicoccum nigrum[END_REF][START_REF] Bensch | Common but different: The expanding realm of Cladosporium[END_REF][START_REF] Setati | Sequence-based Analysis of the Vitis vinifera L. cv Cabernet Sauvignon Grape Must Mycobiome in Three South African Vineyards Employing Distinct Agronomic Systems[END_REF][START_REF] Swett | Characterizing the Brown Spot Pathosystem in Late-Harvest Table Grapes ( Vitis vinifera L.) in the California Central Valley[END_REF][START_REF] Dissanayake | Direct comparison of culture-dependent and culture-independent molecular approaches reveal the diversity of fungal endophytic communities in stems of grapevine (Vitis vinifera)[END_REF], these associations were difficult to interpret due to a lack of knowledge of microbial interactions in natura. Nevertheless, these results show that the combination of network replicates and inference methods permits the identification of apparently robust associations between abundant species, which could be indicative of ecological interactions.

In our study, community-level analyses were found to be more informative, from an ecological perspective, than network-level analyses. We found that the richness, diversity and evenness of fungal communities were significantly higher in organic than conventional vineyards, consistent with the recent findings of [START_REF] Kernaghan | Fungal endophytes of wild and hybrid Vitis leaves and their potential for vineyard biocontrol[END_REF] (but see [START_REF] Castañeda | Effects of agricultural management on phyllosphere fungal diversity in vineyards and the association with adjacent native forests[END_REF]. The cropping system also significantly affected the composition of grapevine foliar fungal communities, as reported in previous studies [START_REF] Castañeda | Effects of agricultural management on phyllosphere fungal diversity in vineyards and the association with adjacent native forests[END_REF][START_REF] Kernaghan | Fungal endophytes of wild and hybrid Vitis leaves and their potential for vineyard biocontrol[END_REF]Pancher et al., 2012;[START_REF] Schmid | Functional and structural microbial diversity in organic and conventional viticulture: Organic farming benefits natural biocontrol agents[END_REF][START_REF] Varanda | Fungal endophytic communities associated to the phyllosphere of grapevine cultivars under different types of management[END_REF]. For instance, Erysiphe necator, the causal agent of grapevine powdery mildew, was significantly more abundant in conventional than in organic plots according to DNA metabarcoding data. These results were consistent with visual assessments of disease symptoms, indicating that, despite their numerous biases, metabarcoding data do contain some quantitative information useful for monitoring plant disease development [START_REF] Jakuschkin | Deciphering the Pathobiome: Intra-and Interkingdom Interactions Involving the Pathogen Erysiphe alphitoides[END_REF]Makiola et al., 2018;[START_REF] Sapkota | Host genotype is an important determinant of the cereal phyllosphere mycobiome[END_REF]. The cause for such contrast in the pathogen abundance is possibly the nature and timing of phytosanitary treatments, but not the dose or number of applications that was similar in the two systems (Table S1). Phytosanitary treatments also influenced several yeast strains, assigned to the genera Vishniacozyma, Sporobolomyces and Filobasidium, that were significantly more abundant in organic plots. These yeast genera are frequently detected on leaf surfaces due to their tolerance of irradiation and they might influence plant growth by producing plant hormone-like metabolites [START_REF] Kemler | Phylloplane yeasts in temperate climates[END_REF]. In addition, Vishniacozyma victoriae (ex Cryptococcus victoriae) was reported as a biocontrol agent of postharvest diseases [START_REF] Lutz | Efficacy and putative mode of action of native and commercial antagonistic yeasts against postharvest pathogens of pear[END_REF]. Other yeasts possessvaluable features of biocontrol agents including killer activities for some Sporobolomyces yeasts [START_REF] Klassen | Antagonistic interactions and killer yeasts[END_REF]. The yeasts Vishniacozyma victoriae and Filobasidium wieringae (ex Cryptococcus wieringae) were also reported as moderate antagonists of several filamentous fungi (Hilber-Bodmer, [START_REF] Hilber-Bodmer | Competition assays and physiological experiments of soil and phyllosphere yeasts identify Candida subhashii as a novel antagonist of filamentous fungi[END_REF]. Future research should investigate the interactions between these yeast species and grapevine foliar pathogens, including powdery mildew.

In the future, we envision that the analysis of microbial interaction networks in the phyllosphere (i.e. the microbial habitat formed by plant leaves [START_REF] Vacher | The Phyllosphere: Microbial Jungle at the Plant-Climate Interface[END_REF][START_REF] Vorholt | Microbial life in the phyllosphere[END_REF])) will serve the prediction of foliar disease risk in crop plants. Plant-associated microbial interaction networks can protect plants against disease [START_REF] Hassani | Microbial interactions within the plant holobiont[END_REF][START_REF] Kemen | Microbe-microbe interactions determine oomycete and fungal host colonization[END_REF]. Resistance to pathogens is mediated by direct antagonistic interactions between the resident microbiota and the invading pathogen species (i.e. the barrier effect; [START_REF] Arnold | Fungal endophytes limit pathogen damage in a tropical tree[END_REF][START_REF] Kamada | Control of pathogens and pathobionts by the gut microbiota[END_REF][START_REF] Kemen | Microbe-microbe interactions determine oomycete and fungal host colonization[END_REF][START_REF] Koch | Socially transmitted gut microbiota protect bumble bees against an intestinal parasite[END_REF][START_REF] Laur | Effectors involved in fungal-fungal interaction lead to a rare phenomenon of hyperbiotrophy in the tritrophic system biocontrol agent-powdery mildew-plant[END_REF] and by indirect interactions due to the activation of the host immune system by the resident microbiota (i.e. the priming effect; [START_REF] Hacquard | Interplay Between Innate Immunity and the Plant Microbiota[END_REF][START_REF] Kamada | Control of pathogens and pathobionts by the gut microbiota[END_REF][START_REF] Perazzolli | Downy mildew resistance induced by Trichoderma harzianum T39 in susceptible grapevines partially mimics transcriptional changes of resistant genotypes[END_REF][START_REF] Ritpitakphong | The microbiome of the leaf surface of Arabidopsis protects against a fungal pathogen[END_REF][START_REF] Vogel | The Arabidopsis leaf transcriptome reveals distinct but also overlapping responses to colonization by phyllosphere commensals and pathogen infection with impact on plant health[END_REF]. The subset of the hostassociated microbial network, consisting of a pathogen and its interacting partners has been termed the pathobiome [START_REF] Brader | Ecology and Genomic Insights into Plant-Pathogenic and Plant-Nonpathogenic Endophytes[END_REF][START_REF] Vayssier-Taussat | Shifting the paradigm from pathogens to pathobiome: new concepts in the light of meta-omics[END_REF]. To better understand and predict disease risk, we should identify the microbial interactions forming pathobiomes [START_REF] Durán | Microbial Community Composition in Take-All Suppressive Soils[END_REF] and the intrinsic network properties that hinder invasion by pathogens [START_REF] Agler | Microbial Hub Taxa Link Host and Abiotic Factors to Plant Microbiome Variation[END_REF][START_REF] Murall | Invasions of Host-Associated Microbiome Networks[END_REF][START_REF] Poudel | Microbiome networks: A systems framework for identifying candidate microbial assemblages for disease management[END_REF]. NGB will require the monitoring in real-time of these properties, based on the automated sequencing on leaf DNA. However, our study shows that statistical network inference, as currently based on a limited sampling effort, generates very few robust hypotheses for microbial interactions, limiting its use to monitoring the disease regulation services provided by the microbiota.

To conclude, here we have demonstrated that microbial networks, automatically inferred from DNA metabarcoding data at the ASV level [START_REF] Callahan | DADA2: High-resolution sample inference from Illumina amplicon data[END_REF] with classical methods of statistical network inference such as SparCC [START_REF] Friedman | Inferring Correlation Networks from Genomic Survey Data[END_REF] or SPIEC-EASI [START_REF] Kurtz | Sparse and Compositionally Robust Inference of Microbial Ecological Networks[END_REF], and then compared using frameworks defined by theoretical ecologists (Pelissier et al., 2018;Poisot et al., 2012), can detect ecosystem change and therefore have a role to play in NGB approaches. Our results suggest that network β-properties were better indicators of change than network α-properties and should be preferred in future developments of NGB. We also showed that keeping the sequence data at the ASV level, rather than aggregating them at higher taxonomic levels, was preferable because it increased the replicability of the networks within a system. In our study, however, inferred networks were highly variable within a system whatever the method of network inference. Network replicates shared more associations than random networks of the same size, but the few shared associations involved only the most abundant ASVs and contained little ecological information on the functioning of the ecosystem. Future research in microbial network inference should therefore improve the replicability and interpretability of networks by, for instance, inferring ecological interaction types rather than positive and negative associations between microorganisms. Mutual information approaches, based on maximal information coefficients (MIC; [START_REF] Reshef | Detecting novel associations in large data sets[END_REF] could overcome this dichotomy although these approaches have not stood out in the inference benchmarkings done to date [START_REF] Hirano | Difficulty in inferring microbial community structure based on co-occurrence network approaches[END_REF][START_REF] Weiss | Correlation detection strategies in microbial data sets vary widely in sensitivity and precision[END_REF]. All functional and ecological knowledge available on microorganisms needs to be gathered in databases [START_REF] Louca | Decoupling function and taxonomy in the global ocean microbiome[END_REF][START_REF] Nguyen | An open annotation tool for parsing fungal community datasets by ecological guild[END_REF][START_REF] Větrovský | GlobalFungi: Global database of fungal records from high-throughput-sequencing metabarcoding studies[END_REF] and integrated into network inference processes. In a study of trophic networks, [START_REF] Bohan | Automated discovery of food webs from ecological data using logic-based machine learning[END_REF] showed that logic-based machine learning is a promising tool to integrate background knowledge to network inference. Future research should investigate the relevance of this approach to microbial network inference. In our study, community-level analyses of DNA metabarcoding data were more statistically powerful than network-level analyses, because many samples were needed to build each network, and this reduced the number of network replicates by comparison with community replicates. The number of samples recommended in the literature for building a single network varies, from 25 [START_REF] Berry | Deciphering microbial interactions and detecting keystone species with co-occurrence networks[END_REF] to 200 [START_REF] Hirano | Difficulty in inferring microbial community structure based on co-occurrence network approaches[END_REF].

Our study shows that networks built from fewer samples (20 in the present case) can nevertheless detect ecosystem change, although we would advise more samples to increase the robustness of the inferred networks. In contrast to network-level properties, all community-level properties detected ecosystem change and provided information important for our understanding of the ecosystem functioning, such as for instance the higher microbial diversity and lower pathogen abundance under organic farming. Community-level analyses should therefore not be discarded in future developments of NGB, that will have to rely on very large DNA metabarcoding datasets combined with functional databases to fully benefit from network-level approaches.

Isabelle Lesur, Adrien Rush, Lucie Zinger and four anonymous reviewers for very helpful comments on the results. We thank the Genotoul sequencing facility (Get-PlaGe) for sequencing and the Genotoul bioinformatics facility (Bioinfo Genotoul) for providing computing and storage resources. We also

Figures and Tables

Figure 1 -Experimental design. Foliar fungal communities were characterized in three conventional (CONV) and three organic (ORGA) vineyard plots by a metabarcoding approach. We analyzed 20 foliar samples per plot. For each plot, we thus obtained 20 community profiles (described in terms of amplicon sequence variants (ASV)) and one association network (inferred either with the SparCC software developed by [START_REF] Friedman | Inferring Correlation Networks from Genomic Survey Data[END_REF] or with the SPIEC-EASI software developed by [START_REF] Kurtz | Sparse and Compositionally Robust Inference of Microbial Ecological Networks[END_REF]. More networks were then obtained by varying network reconstruction parameters (Figure 3). The effects of cropping system (CONV versus ORGA) on the grapevine foliar microbiota were assessed with both community and network α-and β-properties. S4; * p<0.05; **p<0.01; ***p<0.001). (D) Principal coordinate analysis (PCoA) was used to represent dissimilarities in composition between samples, as assessed with the quantitative and (E) binary Jaccard indices. The effect of the cropping system on both βdiversity metrics was significant, as a single effect for the quantitative Jaccard index and in interaction with block for the binary index (Table 2). Green circles, squares and triangles correspond to samples collected in the ORGA1, ORGA2 and ORGA3 plots, respectively. Orange circles, squares and triangles correspond to the CONV1, CONV2 and CONV3 plots, respectively (Figure 1). (F) Log-transformed ratio of ASV relative abundance in CONV plots over that in ORGA plots, for 14 ASVs identified as differentially abundant between cropping systems by DESeq2 analysis followed by Benjamini-Hochberg adjustment [START_REF] Love | Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[END_REF]. 1) were calculated for each network: the number of links (L) and connected components (CC), the network diameter (DIA) and connectance (C) and the mean degree (DEG) and negative link ratio (NLR). The percentage P of ASVs used for network reconstruction had a significant influence on all properties (Table S9), whereas the cropping system did not (Table S7). (C) Principal coordinate analysis (PCoA) represents dissimilarities between networks, measured with the βOS index (Poisot et al, 2012) calculated with the binary Jaccard index. βOS measures the dissimilarity between two networks in terms of the presenceabsence of associations between shared ASVs. The centroids for each cropping system are represented by gray circles. The effect of the cropping system on βOS was significant (Table 4). Networks were inferred with SparCC [START_REF] Friedman | Inferring Correlation Networks from Genomic Survey Data[END_REF]. (A) Associations common to the three network replicates inferred for the organic cropping system (ORGA1, ORGA2, ORGA3) and (B) the three network replicates inferred for the conventional cropping system (CONV1, CONV2, CONV3), regardless of the sign of the association, in the situation in which all ASVs were used for network construction (P=100%). (C) Associations common to the six networks. Networks were inferred with SparCC [START_REF] Friedman | Inferring Correlation Networks from Genomic Survey Data[END_REF]. The number of nodes shared by the network replicates is indicated into brackets. Table 2 -Effect of cropping system -conventional versus organic -on the β-diversity metrics of grapevine foliar fungal communities. Dissimilarities in community composition between samples were assessed with both the quantitative and binary Jaccard indices. The effects of sequencing depth (SD, log-transformed), cropping system (CS) and block (B) on compositional dissimilarities between communities were evaluated using permutational analysis of variance (PERMANOVA), with the number of permutations set to 999.

Dissimilarity index PERMANOVA

Quantitative Jaccard

Variable

thank

  Julie Sappa from Alex Edelman & Associates for English language revision of the first version of the manuscript. We thank the INRA MEM metaprogram (Meta-Omics of Microbial Ecosystems) for financial and scientific support. Sequencing was funded by the INRA MEM MetaBAR project and bioinformatic and statistical analyses were performed as part of the INRA MEM Learn-biocontrol project. Additional funding was received from the LABEX COTE (ANR-10-LABX-45), the LABEX CEBA (ANR-10-LABX-25-01) and INRA EcoServ metaprogram on ecosystem services (IBISC project) and the Aquitaine Region (Athene project, n°2016-1R20301-00007218). CP's PhD grant was funded by the INRA and Bordeaux Sciences Agro (BSA). DBB's grant was funded by the ANR (ANR-17-CE32-0011) and SYNGENTA CROP PROTECTION AG (TK527180). The management of the experimental site was partly funded by the AFB (French Agency for Biodiversity) within the DEPHY network.

Figure 2 -

 2 Figure 2 -Effect of cropping system -conventional (CONV) versus organic (ORGA) -on the αdiversity and β-diversity metrics of grapevine foliar fungal communities. (A) Community richness, defined as the number of ASVs. (B) Community diversity, measured with the inverse Simpson index. (C) Community evenness, measured with Pielou's index. Differences in α-diversity metrics between cropping systems were significant (TableS4; * p<0.05; **p<0.01; ***p<0.001). (D) Principal coordinate analysis (PCoA) was used to represent dissimilarities in composition between samples, as assessed with the quantitative and (E) binary Jaccard indices. The effect of the cropping system on both βdiversity metrics was significant, as a single effect for the quantitative Jaccard index and in interaction with block for the binary index (Table2). Green circles, squares and triangles correspond to samples collected in the ORGA1, ORGA2 and ORGA3 plots, respectively. Orange circles, squares and triangles correspond to the CONV1, CONV2 and CONV3 plots, respectively (Figure1). (F) Log-transformed ratio of ASV relative abundance in CONV plots over that in ORGA plots, for 14 ASVs identified as differentially abundant between cropping systems by DESeq2 analysis followed by Benjamini-Hochberg adjustment[START_REF] Love | Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[END_REF].

Figure 3 -

 3 Figure3-Effect of cropping system -conventional (CONV) versus organic (ORGA) -on the αproperties and β-properties of grapevine foliar fungal networks. (A) Association networks inferred from fungal metabarcoding data with SparCC[START_REF] Friedman | Inferring Correlation Networks from Genomic Survey Data[END_REF]. A total of 60 networks were inferred, corresponding to 2 cropping systems × 3 replicates (blocks) × 10 P values, with P the percentage of most abundant ASVs used for network inference. Only four values of P are shown on the figure. (B) Variations in network α-properties. The following properties (Table1) were calculated for each network: the number of links (L) and connected components (CC), the network diameter (DIA) and connectance (C) and the mean degree (DEG) and negative link ratio (NLR). The percentage P of ASVs used for network reconstruction had a significant influence on all properties (TableS9), whereas the cropping system did not (TableS7). (C) Principal coordinate analysis (PCoA) represents dissimilarities between networks, measured with the βOS index(Poisot et al, 2012) calculated with the binary Jaccard index. βOS measures the dissimilarity between two networks in terms of the presenceabsence of associations between shared ASVs. The centroids for each cropping system are represented by gray circles. The effect of the cropping system on βOS was significant (Table4). Networks were inferred with SparCC[START_REF] Friedman | Inferring Correlation Networks from Genomic Survey Data[END_REF].

Figure 4 -

 4 Figure 4 -Venn diagrams showing the number of fungal associations common to network replicates.(A) Associations common to the three network replicates inferred for the organic cropping system (ORGA1, ORGA2, ORGA3) and (B) the three network replicates inferred for the conventional cropping system (CONV1, CONV2, CONV3), regardless of the sign of the association, in the situation in which all ASVs were used for network construction (P=100%). (C) Associations common to the six networks. Networks were inferred with SparCC[START_REF] Friedman | Inferring Correlation Networks from Genomic Survey Data[END_REF]. The number of nodes shared by the network replicates is indicated into brackets.

Table 1 -List of community-level and network-level α-and β-properties analyzed in the study.

 1 The number of independent observations (N) and the size of corresponding dissimilarity matrices (S) are indicated. The last column indicates if the property varied significantly (Yes/No) with change in the cropping system (CS).

	Property	Definition	Reference	Number of observations	CS
		Community α-properties
	Richness	Total number of amplicons sequence variants (ASVs)	-	N=112	Y
	Diversity (Inverse Simpson)	Effective number of ASVs	Simpson, 1949	N=112	Y
	Evenness (Pielou's J') Evenness in ASV relative abundance	Pielou, 1966	N=112	Y
		Community β-properties
	Compositional dissimilarity (binary Jaccard)	Dissimilarity of composition due to ASV turnover	Jaccard, 1900	S=6216	Y
	Compositional dissimilarity (quantitative Jaccard)	Dissimilarity of composition due to variations in ASV relative abundance	Chao et al, 2006	S=6216	Y
		Network α-properties
	Number of links (L)	Total number of links	-	N=6	N
	Connectance (C)	Fraction of the total number of possible links actually realized	Coleman & Moré, 1983	N=6	N
	Number of connected components (CC)	Number of groups of nodes connected together	Martinez, 1992	N=6	N
	Diameter (DIA)	The longest of all the shortest paths between two nodes	Barabási et al, 2000	N=6	N
	Mean node degree (DEG)	Mean number of links per node	Martinez, 1992	N=6	N
	Proportion of negative links (NLR)	Proportion of links for which the SparCC correlation is negative	Faust et al, 2015	N=6	N
		Network β-properties
	Topological dissimilarity (Schieber's D)	Dissimilarity of global and local network structure	Schieber et al, 2017	S=15	N
	Association dissimilarity (βWN)	Overall dissimilarity of associations	Poisot et al, 2012	S=15	Y
	Association dissimilarity (βOS)	Dissimilarity of associations between shared ASVs	Poisot et al, 2012	S=15	Y
	Association dissimilarity (βST)	Dissimilarity of associations due to ASV turnover	Poisot et al, 2012	S=15	N

Table 3 -Most abundant amplicon sequence variants (ASVs) in grapevine foliar fungal communities according to the cropping system

 3 . The relative abundances (RA, in %) and ranks of ASVs were calculated for all leaf samples (TOTAL; n = 112) and for samples collected from organic (ORGA; n = 55) and conventional plots (CONV; n = 57).

	ASV taxonomic assignment	TOTAL	ORGA	CONV
		Rank	RA	Rank	RA	Rank	RA
	Aureobasidium sp.	1	61.4	1	55.8	1	66.7
	Cladosporium delicatulum	2	6.3	4	6.9	2	5.8
	Filobasidium sp.	3	5.1	2	9.7	9	0.7
	Alternaria sp.	4	4.4	5	3.9	4	5.0
	Epicoccum nigrum	5	4.1	7	2.7	3	5.4
	Cladosporium ramotenellum	6	3.5	3	7	46	<0.1
	Mycosphaerella tassiana	7	3.3	8	1.8	5	4.8
	Didymella sp.	8	1.4	6	2.7	33	0.1
	Erysiphe necator	9	1.1	38	<0.1	6	2
	Vishniacozyma victoriae	10	0.9	9	1.6	17	0.3

Acknowledgements

We thank Lucile Muneret, Andreas Makiola and all members of the ANR NGB Consortium (ANR-17-CE32-0011) for helpful scientific exchanges on next-generation biomonitoring and for their comments on the manuscript. We also thank Gregory Gambetta, Guilherme Martins, Frédéric Barraquand,

Data Accessibility Statement

The raw sequence data were deposited in Dataverse and are available in the FASTQ format at https://doi.org/10.15454/3DPFNJ while the filtered ASV table is available at https://doi.org/10.15454/WOICSE. The code is available as an archive at https://doi.org/10.15454/ZWDFJK.

Author Contributions

CP and DBB performed the bioinformatic and statistical analyses. JV and CV performed the sampling and processed the leaf samples. LD managed the sampling site and provided data on phytosanitary treatments and disease symptoms. JV performed the DNA extractions and amplifications. MB coordinated the work on the fungal mock community and the sequencing of all samples. DB helped interpreting the results in the context of next-generation biomonitoring. CV conceived the study, supervised the analyses and wrote the manuscript with the help of CP and DBB. All authors discussed the results and revised the manuscript.

Table 4 -Effect of cropping system -conventional versus organic -on the β-properties of grapevine foliar fungal networks inferred with SparCC. The D index quantifies the topological dissimilarity between networks [START_REF] Schieber | Quantification of network structural dissimilarities[END_REF] whereas the other three metrics (βWN, βOS and βST), which were calculated with the binary Jaccard index, quantify differences in associations between networks (Poisot et al, 2012). The effect of the percentage P of the most abundant ASVs used for network inference, and the effect of cropping system (CS) on the dissimilarities between networks were evaluated in permutational analysis of variance (PERMANOVA). The number of permutations was set to 999 and permutations were constrained by block.

44

Table 5 -Number of associations shared between network replicates within each cropping systemconventional (CONV) and organic (ORGA) -depending on the method of network inference. Networks were inferred with SparCC [START_REF] Friedman | Inferring Correlation Networks from Genomic Survey Data[END_REF] or SPIEC-EASI [START_REF] Kurtz | Sparse and Compositionally Robust Inference of Microbial Ecological Networks[END_REF], by aggregating or not the ASVs at the genus level, and by including various percentages P of the most abundant ASVs or genera in the network.

The number of shared ASVs or genera between the three network replicates is given into brackets. For every combination of parameters, three random networks having the same number of nodes and links than the three inferred networks were simulated. The pseudo p-value is the probability, estimated with 999 simulations, that the three random networks shared more associations than the three inferred networks (* p<0.05; ** p<0.01; *** p<0.001).

P (%)