Jaime Aguayo

Claude Husson

Emilie Chancerel

Olivier Fabreguettes

Anne Chandelier

Céline Fourrier-Jeandel

Nadine Dupuy

Cyril Dutech

Renaud Ioos

Cécile Robin

Michel Thibaudon

Benoit Marçais

Marie-Laure Desprez-Loustau

Combining permanent aerobiological networks and molecular analyses for large-scale surveillance of forest fungal pathogens: a proof-of-concept

Keywords: Aerobiology, spore dispersal, spore trap, metabarcoding, Hymenoscyphus fraxineus, plant disease surveillance

Forest disease management relies principally on a preventive approach in which epidemiological surveillance plays a crucial role. However, efficient and cost-effective surveillance methods are not currently available for large spatial scales. Yet, aerobiological networks have been set up for several decades in many countries to monitor pollen dispersal and provide real-time assessments of allergenic risk. Here, we suggest that the same approach could be used for the surveillance of forest pathogens. Using molecular methods, we analyzed samples from 12 sites of the French aerobiological network, at different dates. Both metabarcoding by high throughput sequencing (using two different bioinformatics approaches and two markers) and real-time PCR targeting eight important forest pathogens were conducted. To validate the approach, temporal and spatial trends of spore detection were compared with field disease data. The metabarcoding approach demonstrated that many fungal plant pathogens could be found in aerobiological samples. Moreover, five of the eight targeted forest pathogens were detected by real-time PCR, with temporal and spatial trends of spore captures consistent with field data. In particular, Hymenoscyphus fraxineus was detected at high frequency in aerobiological samples in the areas where ash dieback has been present for the longest period of time and at lower frequency in areas with more recent invasion. Spore detection of seasonal pathogens showed a temporal pattern similar to that of disease reports. Overall, our study provides a proof of concept that permanent aerobiological networks combined with molecular methods may provide a useful tool for large-scale surveillance of forest pathogens.

Introduction

Improving disease surveillance for plants, especially at large spatial scales, is a key factor for setting up adequate measures to prevent the establishment and spread of pathogens [START_REF] Parnell | Surveillance to inform control of emerging plant diseases: an epidemiological perspective[END_REF]. In an era of increasing numbers of emerging diseases, more efficient surveillance strategies can foster a quick response to disease outbreaks. In particular, surveillance can be useful to detect pathogens in a given area before they reach a prevalence level that precludes any possibility of eradication or containment [START_REF] Parnell | Surveillance to inform control of emerging plant diseases: an epidemiological perspective[END_REF]. These considerations also apply to forests, which are increasingly threatened by emerging diseases either caused by pathogens of exotic origin, or favored by climate change, or intensification of silvicultural practices [START_REF] Trumbore | Forest health and global change[END_REF]. Moreover, disease control methods are very limited once a pathogen establishes in forests, because fungicides are not economically or environmentally feasible. Minimizing disease risk in forestry thus largely relies on preventive methods whose efficiency depends on accurate risk assessment, for which surveillance data can provide crucial information. For non-native forest pathogens, the spatial extent of ongoing epidemics is of particular importance from both a regulatory and management (including eradication) perspective [START_REF] Filipe | Landscape epidemiology and control of pathogens with cryptic and long-distance dispersal: sudden oak death in northern Californian forests[END_REF]. Large-scale monitoring of major pathogens, including native ones, is also important to provide spatial indicators of forest health, monitor temporal trends, and thus help managers or policy makers implement durable forest management. Various surveillance schemes based on different methods and spatial scales have been proposed or applied in forestry [START_REF] Trumbore | Forest health and global change[END_REF]. It has been for example shown that the accuracy of surveillance data can strongly depend on the underlying sampling methods [START_REF] Hanks | ͘ Reconciling multiple data sources to improve accuracy of large-scale prediction of forest disease incidence[END_REF]. Spatial scale is also a key feature to be considered and often represents the best compromise between disease dynamics and management objectives [START_REF] Meentemeyer | Landscape epidemiology of emerging infectious diseases in natural and human-altered ecosystems[END_REF]. However, the combination of several spatial scales is often required to reconcile accuracy and scale issues [START_REF] Hanks | ͘ Reconciling multiple data sources to improve accuracy of large-scale prediction of forest disease incidence[END_REF]. For example, forest surveillance typically relies on remotesensing methods, often at very large scales [START_REF] Coops | Assessment of Dothistroma needle blight of Pinus radiata using airborne hyperspectral imagery[END_REF], and on specific ground-based surveys at much more local scales. However, producing near-real-time maps of the presence/absence of a disease can be time-consuming and costly, with either remote sensing or classical ground-based methods (surveys carried out by trained observers) [START_REF] Labbé | Pre-existing forests as sources of pathogens? The emergence of Armillaria ostoyae in a recently planted pine forest[END_REF].

Here, we propose a methodology based on aerobiological sampling in a largescale network combined with molecular identification of pathogenic fungi as a potential cost-efficient monitoring scheme contributing to forest health surveillance. To our knowledge, there is currently no perennial dedicated network for the detection of airborne inoculum of pathogenic fungi. However, airborne monitoring networks have been set up for several decades in many countries to monitor allergenic particles (mostly pollens) in the air and provide allergen risk assessments for human populations. In France, for example, the RNSA (Réseau National de Surveillance Aérobiologique) has set up more than 70 volumetric aerobiology traps in urban zones and produces vigilance maps with a grid of roughly 80×80 km, with weekly updates (https://www.pollens.fr/). In our study, the working hypothesis is that large-scale and permanent aerobiology networks used for monitoring allergen particles, such as the RNSA network, can also be used for the surveillance and epidemiology of forest fungal pathogens. This hypothesis is based on the following three premises. First, the same type of particle samplers that are used for pollen monitoring (Hirst-type volumetric air samplers) are routinely used for fungal spore trapping in other settings; only their location is different for pollen monitoring (i.e. placed at several meters high and in urban environments). Second, the traps used for pollen monitoring, thanks to their position at 10 to 20 m above ground level, have a regional representativeness. Indeed, most suspended particles trapped in the air originate from sources located in a radius of several tens of km. For example, in analyses of sources for oak and olive pollen in Spain, the maximum influence of emission sources was observed at 14 km and 37 km from the rooftop pollen trap location, respectively [START_REF] Oteros | Concentric Ring Method for generating pollen maps. Quercus as case study[END_REF]. Finally, currently available molecular methods (targeted or without a priori methods such as real-time PCR or highthroughput sequencing) can provide tools that potentially allow the detection of organisms even at low spore concentrations in air masses. The objective of our study was to demonstrate the feasibility and relevance of using pollenmonitoring networks, combined with molecular biology methods, for the large-scale surveillance of pathogenic fungi, in particular forest pathogens. The first issue was to confirm that forest pathogenic fungi can be detected in aerobiological samples from pollen monitoring networks, despite the location of traps on roof-tops in urban environments (providing regional sampling but with probably very low densities of spores). Although a few studies report the detection of some crop pathogens in rooftop traps [START_REF] Brittain | Detection of Plant Pathogen Spores of Economic Significance on Pollen Trap Slides[END_REF][START_REF] Banchi | DNA metabarcoding uncovers fungal diversity of mixed airborne samples in Italy[END_REF][START_REF] Migliorini | Temporal patterns of airborne Phytophthora spp. in a woody plant nursery area detected using real-time PCR[END_REF], they do not document the presence of forest pathogens. To answer this question, we used two molecular detection methods, a metabarcoding approach and specific real-time PCR. Metabarcoding was implemented with two bioinformatics approaches and two markers targeting all fungal taxa, with a posteriori taxonomical and functional assignments. Specific real-time PCR targeted eight important forest fungal pathogens. The second issue was to assess the epidemiological relevance of detections, i.e. to test whether the frequency of detection of specific pathogens showed temporal or spatial patterns consistent with the occurrence of disease as reported from ground-based surveys.

Here we compared the real-time PCR detections to disease data from the French Forest health service (DSF).

Materials and Methods

Spore sample collection

The RNSA network (French Aerobiology Network) is the French organization in charge of monitoring biological allergen particles (principally pollen, but also fungal spores), and providing real-time information on allergy risks for human populations. The RNSA has set up and maintains a network of Hirst-type volumetric air samplers placed in 70 French cities, distributed all over the country in sites chosen according to climate, botanical and population density criteria. The samplers are located on rooftops, according to a standardized procedure, at a height of 10 to 20 m to take air samples representative of the local air masses. In this study, we used 12 Hirst-type volumetric air samplers (located in 12 different French urban areas) (Table 1). In order to maximize the chances of detection, sites were selected according to the current distribution of forest diseases among those most frequently reported in mainland France and caused by fungi with aerial dispersal. Air samples were collected during the year 2016 and sampling dates for analysis were selected for each pathogen by focusing on the most likely period of spore dispersal according to the literature (Table 1).

Sample processing

A total of 154 samples were processed (Table 1). The samples consisted of microscope slides prepared by the RNSA for pollen counts, according to a standardized procedure (European standard: CEN/TC 264/WG 39). Briefly, Hirsttype air samplers are active suction traps, functioning with a vacuum pump providing a constant input airflow of 10 l per min. Upon impact, airborne particles adhere to a Vaseline-coated tape (Melinex tape) mounted on a drum that rotates at a constant speed (2 mm per hour). Recordings are generally run continuously for seven days. At the end of each recording period, the tape is divided into sections, each section representing one day (a 24 h sample corresponds to 48 mm). Each tape section is mounted on a microscope slide stained with fuchsine dissolved in lactic acid and covered with a glass cover slip. These microscope slides, each corresponding to a 24 h period and 14.4 m3 of air, were the starting material for this study.

The procedure for DNA extraction, selected after a preliminary experiment comparing several laboratory treatments (data not shown), was performed as follows. Mounted microscope slides were placed flat on a constant heater set at 65°C for 4 to 5 min, to remove the glass cover slip. The sticky tape was recovered with forceps and cut into small pieces with sterile scissors. The tape pieces were then transferred into a 1.5 ml tube containing Lysing Matrix C and 1 mm silica spheres (MP Biomedicals, Illkirch-Graffenstaden, France), filled with 400 µl of AP1 buffer and 4 µl of RNase A (both from the DNeasy plant mini kit Qiagen, Hilden, Germany).

The content of the tube was ground twice for 60 s using a high-speed FastPrep®24 (MP Biomedicals) homogenizer set at 6 shakes s -1 . Total genomic DNA was then extracted from the mash (tape + lysis buffer) with the Qiagen DNeasy plant mini kit, following the manufacturer's specifications. Total DNA was eluted in 100 µl AE buffer (Qiagen) and kept at -20°C until use as template for real-time PCR and/or metabarcoding analysis. This protocol was strictly applied in the three laboratories where DNA extractions were performed (please refer to Supplementary Material S1).

High-throughput sequencing metabarcoding

Metabarcoding was performed by sequencing the internal transcribed spacer (ITS) of the fungal ribosomal DNA using an Illumina MiSeq Platform. A subsample of 96 microscope slides taken from the different locations and dates was used for the high-throughput sequencing (HTS) metabarcoding study (Table 1). Two primer pairs, Targeting the two spacers of the ITS region, ITS1 and ITS2, were used. Both loci (ITS1 and ITS2) were sequenced because results of fungal community composition may differ significantly according to the genomic region targeted [START_REF] Aguayo | Assessment of Passive Traps Combined with High-Throughput Sequencing To Study Airborne Fungal Communities[END_REF]. The ITS1 region was amplified using the primer pair ITS1F / ITS2 [START_REF] Gardes | ITS primers with enhanced specificity for basidiomycetes -application to the identification of mycorrhizae and rusts[END_REF][START_REF] White | Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics[END_REF], and the ITS2 region was amplified using the primer pair ITS86 / ITS4 [START_REF] White | Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics[END_REF][START_REF] Turenne | Rapid Identification of Fungi by Using the ITS2 Genetic Region and an Automated Fluorescent Capillary Electrophoresis System[END_REF]. Mix and PCR conditions are presented in Supplementary data S2. Protocols used for HTS library preparation and sequencing are detailed in Supplementary Material S3.

Bioinformatics analysis of metabarcoding data

Bioinformatics analyses were performed using two alternative approaches, here named OTU (Operational Taxonomic Units) and ASV (Amplicon Sequence Variants). In an OTU approach, sequences are clustered together based on a fixed similarity threshold (in fungi, usually 97%), that determines their similarity to each other [START_REF] Mysara | Reconciliation between operational taxonomic units and species boundaries[END_REF]. In contrast, an ASV approach does not cluster sequences using a defined similarity threshold, but uses all sequences after correcting sequencing errors. This approach is based on the hypothesis that variations in sequences due to amplification and sequencing errors can be discriminated from true biological variants using bioinformatics analytical methods; the ASV approach is thus expected to provide high resolution (down to the level of single-nucleotide differences) and reproducible results [START_REF] Callahan | DADA2: High-resolution sample inference from Illumina amplicon data[END_REF]. However, there is still no consensus on which method should be preferred [START_REF] Pauvert | Bioinformatics matters: The accuracy of plant and soil fungal community data is highly dependent on the metabarcoding pipeline[END_REF], so we tested and compared both.

For analyses following the OTU approach, sequences were analyzed using PIPITS [START_REF] Gweon | PIPITS: an automated pipeline for analyses of fungal internal transcribed spacer sequences from the Illumina sequencing platform[END_REF], an automated pipeline adapted to fungal ITS sequences from Illumina MiSeq platforms that includes several third-party tools (https://github.com/hsgweon/pipits). For the ASV approach, reads were analyzed with the ITS pipeline workflow (https://benjjneb.github.io/dada2/ITS_workflow.html) implemented in the DADA2 (version 1.8.0) program [START_REF] Callahan | DADA2: High-resolution sample inference from Illumina amplicon data[END_REF]. For both OTUs and ASVs, taxonomic assignment was performed using the UNITE reference database [START_REF] Abarenkov | The UNITE database for molecular identification of fungi -recent updates and future perspectives[END_REF], using the default algorithms proposed in PIPITS and DADA2. An additional manual OTU table filtering was performed using the BLAST algorithm which was run against the UNITE reference database. This was performed to check and correct, if necessary, OTU and ASV assignments at taxonomic levels higher than species. Detailed bioinformatics steps are presented in the Supplementary Material S4.

Only OTUs and ASVs with more than 10 sequences were conserved for downstream analysis. Four datasets were consequently constructed (two markers × two bioinformatics approaches). Because the metabarcoding approach is used here as a surveillance tool, we pooled OTUs or ASVs with the same taxonomical assignment (for example, several OTUs assigned to Alternaria sp. were considered as one OTU in the respective datasets).

FUNGuild [START_REF] Nguyen | FUNGuild: An open annotation tool for parsing fungal community datasets by ecological guild[END_REF] was used to assign functional guilds to datasets generated from both approaches. Fungal taxa were considered plant pathogens whenever the guild "plant pathogen" was mentioned, either as principal or secondary functional guild͘ . This information was added to the four datasets previously created.

Analysis of metabarcoding datasets

The four datasets (two markers × two bioinformatics approaches) were compared for their ability to detect the eight forest pathogens selected in this study (Table 1). The datasets were also compared for their ability to detect the fungal pathogens that are most frequently reported in the DSF database, i.e. a list of 10 taxa representing approximately 75% of all forest disease reports (Table 3). To confirm the identity of these pathogens, the final validation of the taxonomic assignment of the corresponding sequences of OTUs or ASVs (after automatic taxonomic assignments using the UNITE database, as described previously) was performed manually using the BLAST nucleotide sequence tool, with a similarity threshold of 97% with the NCBI database (GenBank).

Specific real-time PCR

Real-time PCR (qPCR) with specific primers were used to detect eight taxa (six species and two species complexes). They corresponded to important forest pathogens (either native or introduced) in France with known or putative aerial dispersal: H. fraxineus, M. larici-populina, C. parasitica, Dothistroma pini and D. septosporum, E. alphitoides, the Heterobasidion annosum complex (comprising H. abietinum, H. annosum s.s., H. irregulare, H. occidentale and H. parviporum), and Armillaria spp. (Table 2). For each targeted taxon, a subset of trap locations and dates were used (Table 1). Detailed information on PCR conditions and parameters for each targeted taxon are presented in Supplementary Material S2. qPCR results were analyzed qualitatively, i.e. as presence/absence of the target fungal species. No Ct threshold (i.e. limit of sensitivity of the assays) was fixed as [START_REF] Grosdidier | Detection of plant pathogens using real-time PCR: how reliable are late Ct values?[END_REF] have shown that late-cycle threshold values (Ct) can be reliably rated as positive, especially when the target DNA is present in low quantities. An estimation of spore concentrations in positive samples was nevertheless computed for H. fraxineus, for which aerobiological samples were compared with spore dilutions in qPCR, as described by [START_REF] Grosdidier | Tracking the invasion: dispersal of Hymenoscyphus fraxineus airborne inoculum at different scales[END_REF].

Comparisons of spore detection and disease data

As explained before, the selection of sites and dates was done in order to maximize detection of the targeted fungi (i.e. spatial and temporal samplings within the known spatial range and spore emission period of these fungi). However, in order to test the epidemiological significance of detection data, some spatial and temporal variations were included. For each fungus, samplings were analyzed at different dates extending beyond the expected spore emission period. In addition, several sites, potentially varying in disease severity, were analyzed. For example for Hymenoscyphus fraxineus, the agent of ash dieback, which was first reported in North-East France in 2008 and is still spreading southwest, six sites were analyzed, two in the area with earliest presence and highest severity of disease, two in recently affected areas, and two outside beyond the invasive front. Disease data were extracted from the French National Forest Health Department database (Département de la santé des forêts, DSF), initiated in 1989 and containing more than 50,000 reports. Logistic regressions performed with the R environment (version 3.5.1) were used to assess (i) the relationship between the detection frequency in the different sites obtained by qPCR (during the suitable period for leaf pathogens) and local disease indices (number of years of disease presence for H. fraxineus and computed report rate in the DSF database for all other pathogens) and (ii) the relationship between monthly spore detection by qPCR and monthly disease reports in the DSF database.

Results

Metabarcoding

Raw sequence reads have been deposited in the NCBI Sequence Read Archive (SRA) under accession number PRJNA607907 and PRJNA608245 for ITS1 and ITS2, respectively. A summary of results from the metabarcoding analysis are presented in Table 4. Similar proportions of OTUs and ASVs, around 21-26%, were assigned as potential plant pathogens (principal or secondary functional guild) based on FUNGuild, with both bioinformatics methods and the two markers. Analyses based on the OTU approach resulted in 1,541,756 (34.35% of the total number of reads) and 2,469,881 (50.44%) valid reads assigned as potential plant pathogens, for ITS1 and ITS2, respectively. They corresponded to 225 and 267 OTUS 21.63% and 22.98% of the total number of OTUs) after pooling for ITS1 and ITS2, respectively. In the ASV approach, there were 1,914,626 (49.57%) and 2,297,207 (52.63%) reads assigned to potential plant pathogens, corresponding to 293 (25.15% of the total number) and 371 (25.46%) ASVs for ITS1 and ITS2, respectively. Detailed OTU and ASV tables are presented in the Supplementary material S5.

In terms of taxonomic assignments, differences were observed between the four datasets, because only 15.2% (369 pooled OTUs and ASVs) were common to all datasets and only 9% (57 pooled OTUs or ASVs) when considering only the plant pathogens. For instance, Hymenoscyphus fraxineus and Melampsora laricipopulina were detected only with the ITS2 marker using either approach (OTU or ASV) whereas Armillaria spp. was detected only with ITS1 using either approach (Table 3). In other cases, results differed according to the bioinformatics approach. For example, Nothophaecryptopus gaeumannii was only detected using the OTU approach, with the ITS1 and ITS2 regions. Discrepancies between bioinformatics methods and markers were especially frequent at the species level.

For example, Armillaria spp. was assigned to A. socialis (i.e. syn. A. tabescens) and A. mellea in the OTU approach, whereas it was associated with A. tabescens following the ASV approach (in addition to assignments at the genus level). Within the Heterobasidion annosum complex, sequences were assigned to H. abietinum in the OTU approach (with both markers), but to H. occidentale and H. orientale, an American and Asian species [START_REF] Chen | Phylogeny, divergence time estimation, and biogeography of the genus Heterobasidion (Basidiomycota, Russulales)[END_REF] respectively, in the ASV approach. The manual validation of taxonomic assignments at the species level confirmed the detection of Erysiphe quercicola, H. fraxineus, M. larici-populina, Gymnopus fusipes, and N. gaeumannii. Other sequences, with an automatic assignment (with the default algorithms of PIPITS and DADA2) to various species of Armillaria, Heterobasidion, Erysiphe, Hymenoscyphus and Melamspora, could not be assigned to the species level since the corresponding sequences matched several species at the 97% threshold level in the BLAST analysis. For example, the sequences attributed to H. occidentale in the OTU or ASV approach had a similarity >97% not only with that species but also with H. annosum or H. parviporum. These results are reported in detail in Supplementary Material S6.

As a whole, among the dozen of most important forest pathogenic fungi in France, two genera (Armillaria and Heterobasidion) and five species (Erysiphe quercicola, Gymnopus fusipes, H. fraxineus, M. larici-populina, N. gaeumannii) were detected (Table 3).

qPCR for targeted forest pathogens

Five of the eight targeted taxa were detected: Armillaria spp., the H. annosum complex, E. alphitoides, H. fraxineus and M. larici-populina. Armillaria spp. was detected in the three analyzed sites, with the highest frequency in late summer (August-September) (Figure 1). Highest frequencies were detected in Mont-de-Marsan, but also in Clermont-Ferrand, in the Massif Central, both areas with high disease report rates and very low at Bordeaux, more distant from the Landes forest (Figure 2a, likelihood Chisq p < 0.001). Positive cycle thresholds (Ct) were comprised between 32 and 36 (34.61 ± 1.73, mean ± standard deviation (SD)). E. alphitoides, the causal agent of oak powdery mildew, was detected in six of the seven selected study sites. However, its detection frequency was not related to the disease report density (Figure 2a, likelihood Chisq p= 0.894). The first detections of oak powdery mildew spores occurred in May, with a peak in June-July-August then a decrease in September. The lowest Ct values (between 34 and 35, whereas the average for positive detections was 37.17 ± 0.86; mean ± SD) were observed in late June-early July in Nevers and Nancy. The monthly detection frequencies followed the seasonal trend of disease reports (Figure 1, likelihood Chisq p <0.001). Positive detections of H. fraxineus, the causal agent of ash dieback, showed a peak during the early summer (June-July), and then decreased in August (Figure 1). Samples from all the other dates were negative. Positive detections were in good agreement with the known distribution of the disease at the time of sampling (likelihood Chisq p< 0.001): 100% of detection for the two spore traps located in the region where the disease had been reported for the longest time and where the regional host density is the highest (traps located in Nancy and Besançon) (Figure 2b). Lower detection rates were observed in areas only recently affected by the disease and there was no detection beyond the disease front in Bordeaux and Mont-de-Marsan (Figure 2b). Ct values ranged from 37 to 38 for Nancy and Besançon and from 39 to 40 for the other positive sites (mean ± SD: 39.15 ± 1.70), corresponding to approximately 100 and 10 spores/day (7 and 0.7/m 3 air), respectively, based on qPCR efficiency and positive controls. The H. annosum complex was detected throughout the whole sampling period at the three tested sites (Figure 2b) with a frequency related to disease report rate (likelihood Chisq p < 0.001), with peaks in April, June and August (mean Ct values ± SD: 35.73 ± 1.1) (Figure 1). M. larici-populina, the agent of poplar rust, was detected at frequencies ranging from 35% to 100% in four out of the five study sites in northeastern and central France where the disease has been commonly reported (Figure 2b, relation with disease report rate not significant, likelihood Chisq p = 0.466). First positive detections of M. larici-populina occurred in May, with increasing frequencies until August, then decreasing in September (mean Ct values ± SD: 29.72 ± 1.99. This seasonal pattern was consistent with disease reports (Figure 1, likelihood Chisq p < 0.001). Cryphonectria parasitica and the two species of Dothistroma were never detected.

Discussion

Our study provides a proof of concept that the use of permanent airborne particle monitoring networks (such as those used for pollen monitoring) combined with molecular methods could be an efficient tool for large-scale monitoring (not yet speaking of early warning) of at least some forest pathogens. More precisely, we showed that: (i) DNA extraction and amplification can be carried out on slides already prepared for microscopy making it possible to perform the pathogen analysis in addition to the pollen analysis; (ii) many pathogenic fungi can be detected on these slides as shown by metabarcoding; (iii) qPCR on DNA obtained from slides was sufficiently sensitive to detect five out of eight targeted taxa among the most frequent European forest pathogens; and (iv) the spatial and temporal patterns of detection by qPCR in aerobiological samples were in good agreement with spatiotemporal disease data.

Although spore traps have been largely used in epidemiology, so far, most studies have been performed at local scales (i.e. with spore traps usually at ground level, within or at proximity of forests stands or crop fields). The few studies using air samplers set up on rooftops or towers for the detection of plant pathogens, were generally performed at few sites and dates, and none reported the detection of forest pathogens [START_REF] Brittain | Detection of Plant Pathogen Spores of Economic Significance on Pollen Trap Slides[END_REF][START_REF] Banchi | DNA metabarcoding uncovers fungal diversity of mixed airborne samples in Italy[END_REF][START_REF] Migliorini | Temporal patterns of airborne Phytophthora spp. in a woody plant nursery area detected using real-time PCR[END_REF]. Our study takes a step further by clearly positioning this approach in an epidemiological surveillance perspective using an established network of traps, such as those operated for pollen monitoring, and by providing data demonstrating its feasibility and applicability for forest pathogens. Applications may also be found for other airborne plant pathogens, since crop and fruit tree pathogens were detected as well (Supplementary table S5). However, our study also identified several pitfalls and limitations of the proposed methodology, especially for the application of metabarcoding in epidemiology and surveillance. We discuss these points more thoroughly below.

Among the eight taxa targeted in our study, five were detected. In contrast, Dothistroma septosporum, D. pini and Cryphonectria parasitica were never detected, either with HTS or with qPCR. For Dothistroma spp., a first explanation for the lack of detection by qPCR could be an insufficient sensitivity of the method, because only primers targeting single-copy genes, designed for detection in infected plant tissue and not environmental samples, were available (Ioos et al., 2009a). Another possible explanation includes unfavorable weather conditions for the disease in 2016, limiting production and release of aerial inoculum. However, the mode of dispersal of these species may be the critical factor. Both Dothistroma spp. and Cryphonectria are deemed to disperse their spores mostly at short distances by rain splash, in the form of conidia emitted in slimy masses [START_REF] Mullett | A reconsideration of natural dispersal distances of the pine pathogen Dothistroma septosporum[END_REF][START_REF] Rigling | Cryphonectria parasitica, the causal agent of chestnut blight: invasion history, population biology and disease control[END_REF]. Spore concentrations in the air decrease rapidly with distance from the source [START_REF] Aylor | Spread of plant disease on a continental scale: role of aerial dispersal of pathogens[END_REF]. Thus, only fungal species able of long distance aerial dispersal are expected to be found in aerobiological trappings such as those used here. Generally, spores are small particles that are theoretically easily transported even by weak winds [START_REF] Wilkinson | Modelling the effect of size on the aerial dispersal of microorganisms[END_REF]. However, transport of fungal spores in air masses first requires carrying the spores through the boundary layer [START_REF] Pringle | Reaching the wind: Boundary layer escape as a constraint on ascomycete spore dispersal[END_REF]. Fungal species vary enormously in their ability to be released from fruiting bodies or sporophores and propelled into the air, which may explain differences in their dispersal abilities. In particular, active mechanisms for basidiospores and ascospores [START_REF] Pringle | Reaching the wind: Boundary layer escape as a constraint on ascomycete spore dispersal[END_REF] are generally more effective than in conidia, explaining, in part, much higher concentrations of the former than the latter in the air [START_REF] Lacey | The aerobiology of conidial fungi[END_REF]. In the case of C. parasitica, the timing of ascospore production is not precisely known [START_REF] Rigling | Cryphonectria parasitica, the causal agent of chestnut blight: invasion history, population biology and disease control[END_REF], thus our sampling may have missed those periods. Conversely, several species that release only or preferentially ascospores or basidiospores were successfully detected in our study, using either qPCR or metabarcoding. This was especially true for H. fraxineus, for which analyzed slides were mainly selected during the known period of ascospore production in summer [START_REF] Grosdidier | Tracking the invasion: dispersal of Hymenoscyphus fraxineus airborne inoculum at different scales[END_REF] and N. gauemannii, the agent of Douglas-fir Swiss needle cast (Stone et al., 2008). The relatively high detection rates of Heterobasidion spp. and Armillaria spp. may also be explained by their dispersal via basidiospores. For both genera however, the PCR primers were not species-specific, hence, the positive signal may be related to several species in each complex, potentially affecting a wide range of tree species. Using the metabarcoding approach, we also detected several other Basidiomycota fungi that are pathogenic to trees (and likely dispersed by basidiospores) such as Gymnopus fusipes or Ganoderma sp. The two other successfully detected taxa by both qPCR and metabarcoding in our study belong to the powdery mildew and rust fungi, E. alphitoides/E. quercicola and Melampsora larici-populina.In the case of the poplar rust disease, its regional dispersal via urediniospores (asexual spores) has been demonstrated elsewhere [START_REF] Xhaard | Disentangling the genetic origins of a plant pathogen during disease spread using an original molecular epidemiology approach[END_REF]. Other species in these groups were observed in the metabarcoding datasets such as M. populnea and M. epitea, aspen and willow rusts, Puccinia graminis, P. hordei and P. coronata, cereal rusts, Pucciniastrum aerolatum, spruce cone rusts, Erysiphe platani, plane powdery mildew (see Supplementary Table S5). Although dispersal during the growing season is mainly associated with conidia, the long-distance dispersal of powdery mildews as well as rust fungi is well documented [START_REF] Hovmøller | Clonality and long-distance migration of Puccinia striiformis f.sp. tritici in north-west Europe[END_REF].

The main value of metabarcoding is its theoretical comprehensiveness for a given set of primers. Here, we used primers targeting the ITS genomic region of fungi, widely used in this kind of studies [START_REF] Aguayo | Assessment of Passive Traps Combined with High-Throughput Sequencing To Study Airborne Fungal Communities[END_REF][START_REF] Banchi | DNA metabarcoding uncovers fungal diversity of mixed airborne samples in Italy[END_REF]. With this approach, we could demonstrate that many plant, especially forest, pathogens with aerial dispersal can be potentially detected in traps located on rooftops in urban environments. Some examples in our study are N. gaeumannii, Lophodermium spp., agents of pine needle cast, M. populnea and M. epitea, Erysiphe arcuata and Phyllactinia fraxini, agents of powdery mildew on hornbeam and ash, respectively, Fistulina hepatica, a wood-decaying fungus of hardwood trees and Mycosphaerella punctiformis, a foliar endophyte of oaks. Another example is Diplodia seriata, the agent of black dead arm of grapevine, with corresponding sequences frequently detected from the trap located near Bordeaux, close to vineyards. However, our study confirmed the increasingly recognized limitations of metabarcoding approaches, especially when the technique is used for surveillance purposes, due to the occurrence of false negatives and false positives, in part explained by biases and artifacts linked to the markers and bioinformatics approaches [START_REF] Aguayo | Assessment of Passive Traps Combined with High-Throughput Sequencing To Study Airborne Fungal Communities[END_REF][START_REF] Pauvert | Bioinformatics matters: The accuracy of plant and soil fungal community data is highly dependent on the metabarcoding pipeline[END_REF][START_REF] Piper | Prospects and challenges of implementing DNA metabarcoding for high-throughput insect surveillance[END_REF]. For example, no Armillaria spp. could be detected with primers targeting ITS2, or M. larici-populina and Hymenoscyphus with ITS1. Furthermore, only a small proportion of the detected taxa were common to the two used bioinformatics approaches.

The most problematic issue with metabarcoding at its current stage of development from an epidemiological perspective lies in the inaccurate taxonomic assignment at the species level. The problem is crucial when dealing with emerging or regulated organisms. Our study confirmed that results obtained by computational classification methods should be confirmed by manual analyses. For example, our final metabarcoding datasets, especially with the ASV method and ITS2, included several species not known to occur in Europe (Heterobasidion occidentale, H. orientale, Hymenoscyphus dehlii, Melampsora occidentalis), including species on quarantine lists, which are carefully monitored and thus unexpected to be found in environmental samples (i.e. F. circinatum and T. indica). Detailed analysis of sequences using the BLAST algorithm strongly suggested that these results are very likely a result of incorrect species assignments (Supplementary material S6). Despite these limitations, metabarcoding, at its present state of development, may be useful in plant pathogen surveillance if its results (i.e. for emerging and/or quarantine pathogens) are systematically confirmed by duly validated species-specific PCR assays.

With qPCR, we were able to detect H. fraxineus at very low concentrations of less than one spore per cubic meter of air on average. However, coupling high sensitivity and high specificity may be problematic in some practical cases. The ITS region allows high sensitivity due to the presence of multiple copies per genome, but may be less polymorphic than single copy genes, that better discriminate closely related species [START_REF] Feau | Finding single copy genes out of sequenced genomes for multilocus phylogenetics in non-model fungi[END_REF]. As example, the primers used for specific detection of E. alphitoides may amplify DNA from powdery mildew species in the same clade, sharing highly similar ITS sequences but not affecting oaks [START_REF] Chen | Phylogeny, divergence time estimation, and biogeography of the genus Heterobasidion (Basidiomycota, Russulales)[END_REF]. On the other hand, primers targeting single copy genes to differentiate species within the Heterobasidion annosum complex, designed for diagnostic purposes [START_REF] Ioos | ͘ Multiplex real-time PCR assays for the detection and identification of Heterobasidion species attacking conifers in Europe[END_REF], failed to detect any target in aerobiological samples (results not shown).

Epidemiological monitoring for surveillance purposes should be carried out at the same spatial scale as the ultimate control measures. In crops, disease control still depends strongly on fungicide applications; therefore, farmer-led monitoring at the field level appears to be an adequate spatial level [START_REF] Carisse | Grape powdery mildew (Erysiphe necator) risk assessment based on airborne conidium concentration[END_REF]. This strongly differs from disease control strategies in forests, which largely rely on a preventive approach. Therefore, risk analysis, especially hazard mapping, at a regional scale is a more relevant approach [START_REF] Meentemeyer | Landscape epidemiology of emerging infectious diseases in natural and human-altered ecosystems[END_REF]. Our study is a first step toward this general strategy. Our results showed that real-time PCR detection in aerobiological samples is consistent with field disease data for several forest pathogens, in particular temporal patterns for the seasonal diseases (oak powdery mildew and poplar rust), and spatial patterns for the chronic diseases (root and butt rot). Interestingly, temporal patterns of spore detection for seasonal diseases seemed to precede the observation of damage. The lack of a significant relationship between the detection frequency (for year 2016) and spatial disease data from DSF for seasonal diseases could be explained by a year effect, since data from the DSF database are cumulated for the 1989-2016 period. Moreover, for oak mildew, both the host and the pathogen are abundant throughout France, not only in the forests surveyed by the DSF, but also in hedges or landscape trees. The example of H. fraxineus is especially interesting because our findings strongly suggest that aerobiology methods can be used to study the spatial spread of invasive alien pathogens (with natural aerial dispersal). This approach may provide not only qualitative results (presence/absence), but also quantitative assessments, as permitted by qPCR assays. Our preliminary results showed that early Ct values for M. larici-populina, E. alphitoides and were obtained in situations of high-expected spore load from disease data.

From a practical point of view, pollen-monitoring networks, present in many countries, offer a remarkable permanent, well maintained, standardized, continuously operating system, providing data on a weekly basis. Although these systems rely on traps located in urban areas due to their primary focus on monitoring human allergens, our study showed that several forest and other plant pathogens could also be detected. These findings are in agreement with theoretical considerations and empirical data on the regional representativeness of rooftop samplers, with the main contributing sources of pollen and spores occurring in a range of a few tens of kilometers around the trapping site [START_REF] Oteros | Concentric Ring Method for generating pollen maps. Quercus as case study[END_REF].

The same slides can thus also be used for both pollen and fungal spores, whatever their detrimental effect as human allergens, plant, forest or animal pathogens. Because we were able to extract and amplify fungal DNA even after microscopic analysis, using past trap samples may even allow the study of past temporal trends on preserved slides.

Here, we demonstrated the usefulness of a network of rooftop air samplers, located 50 to 100 km apart, for large-scale monitoring of some forest pathogens with aerial dispersal. Such networks already exist for pollen allergy risk assessment, but pollen counts still mostly rely on the traditional microscopic observation method, as performed in the French RNSA network. Our study confirms the high potential of molecular techniques, combined with pollen traps, to characterize airborne organisms, as a less labor-intensive and more accurate method than human observation. In particular, qPCR showed high sensitivity and may make it possible to draw real-time maps of pathogen distributions, including some quantification, without ground surveys. At the current stage of implementation, the metabarcoding approach still has limitations in taxonomic assignment precluding a direct epidemiological application (especially for early warning) but rapid improvements in long reads sequencing and development of bioinformatics methods may offer interesting prospects [START_REF] Zahariev | Cluster oligonucleotide signatures for rapid identification by sequencing[END_REF], Tedersoo & Anslan, 2019). Further improvements may come from a better understanding of fungal spore dispersal (from spore discharge to spore deposition, including atmospheric transport, Redondo et al 2020) and inoculum-disease relationships, which can help design efficient surveillance schemes combining different methods at different spatial scales.

Figure 2a .

 2a Figure 2a. Results for Armillaria sp., Cryphonectria parasitica, Dothistroma sp. (D. pini and D. septosporum), and Erysiphe alphitoides. For E. alphitoides frequencies (pie charts) were computed only during the dispersal period (May-August). Maps

Figure 2b .

 2b Figure 2b. Results for Heterobasidion sp. (H. annosum complex), Hymenoscyphus fraxineus and Melampsora larici-populina. Frequencies (pie charts) were computed only during the dispersal period (June-August for H. fraxineus, May-September for M. larici-populina). Maps show, for each pathogen in the DSF database the density of reports during the 1989-2016 period, except for the H. fraxineus map, which shows the progression of the first reports of ash dieback (dark green in 2008-2010, medium green in 2011-13, light green in 2014-15, white, not yet present in 2016). DSF surveys are not systematically performed each year; we thus pooled the data from 1989 to 2016, because these data are accurate representations of the regions with the greatest frequency of pathogen occurrence.

Table 1 .

 1 Study sites, analyzed sampling dates (in 2016) and target species. N = number of samples (1 slide per sampling date) analyzed using qPCR for pathogens listed in the target species column or by metabarcoding (HTS). Asp, Hsp, Cp, Dot, Ea, Hf, Mlp stand for Armillaria spp., Heterobasidion spp., Cryphonectria parasitica, Dothistroma pini and D. septosporum, Erysiphe alphitoides, Hymenoscyphus fraxineus, and Melampsora larici-populina, respectively. The total number of slides analyzed was 154 for real-time PCR (qPCR) assays and 96 for high-throughput sequencing (HTS) metabarcoding.

	Sites	Latitude Longitude Dates	Target	N
	Besançon	47°14'N 6°12'E	21, 22, 26,27 June; 1-6 July	Hf
			22, 26 June; 1, 3, 5 July	HTS	5
	Bordeaux	44°48'N 0°74'W	8, 23 from Feb. to Dec.	Asp, Hsp, Cp, Ea, Hf, Mlp
			8, 23 from Feb. to Sept.	HTS
	Castres	43°36'N 2°15'E	1, 8, 26 May; 4, 11 June	Cp	5
			8, 26 May; 5, 11 June	HTS	3
	Clermont-	45°47'N 3°05'E	16 Feb.; 1, 16 from Mar. to Sept.	Asp, Hsp, Ea
	Ferrand		16 Feb.; 1, 16 from Mar. to Sept.	HTS
	Mont de	43°89'N 0°30'W
				Hf	8
			1, 16 from Mar. to Aug; 1 Sept.	HTS
	Nancy	48°42'N 6°12'E	21, 22, 26, 27 June; 1-6 July	Ea, Hf , Mlp
			21, 26 June; 1, 3, 5 July	HTS	5
	Nevers	47°00'N 3°09'E	13, 29 May; 5-9, 17, 24, 30 June	Dot, Ea, Hf, Mlp
			6 Sept.
			13 May, 5, 8, 17, 30 June; 6 July	HTS	6
	Orléans	47°54'N 1°54'E	1-8 July, 23-31 Aug.	Cp, Ea, Hf, Mlp
			1, 3, 5, 7 July	HTS	9
			23, 25, 27, 29, 31 Aug.
	Pau	43°18'N 0°62'W	13, 14, 29, 30 May; 5-9, 24, 25 June;	Dot, Cp
			17, 30 July; 6, 7 Sept.
			14, 29 May; 7, 8, 24 June; 30 July	HTS	6
	Périgueux	45°12'N 0°44'E	1, 2, 14, 15 May ; 4, 5, 11, 12, 26, 27	Cp
			June
			15 May; 4, 11, 12, 26 June	HTS	5
	Rennes	48°06'N 1°80'W	12-14, 28-30 May; 4, 5, 8, 9, 23, 25	Ea, Cp
			June, 16, 17, 29, 30 July; 5-7 Sept.
			12,14, 30 May, 4, 9, 24 June	HTS	9
			16, 17 July; 5 Sept.
	Troyes	48°18'N 4°30'E	18-28 Aug.	Mlp, Ea
			18, 20, 23, 26 Aug.	HTS	4

1, 16 from Mar. to Aug.; 1 Aug. Asp, Hsp, Dot Marsan 16 May; 1, 16 June; 1, 16 Jul.; 1-16 Aug.; 1 Sept.

Table 4 .

 4 General statistics for each bioinformatics approach/marker (operational taxonomic units [OTU] or amplicon sequence variants [ASV])/internal transcribed spacer, ITS1 or ITS2)

		OTU ITS1	OTU ITS2	ASV ITS1	ASV ITS2
	Total number of valid reads	4 488 286	4 896 408	3 862 486	4 365 026
	Number of reads
	Range per sample Mean ± SD	131 -73 419 46 751 ± 12 706	7016 -88 858 51 004 ± 7 559	103 -69 284 40 234 ± 11 237	6 633 -81 304 45 469 ± 6850
	Number of taxa (OTU/ASV)
	Range per sample Mean ± SD	11 -299 92 ± 53	12 -318 101 ± 56	4 -282 92 ± 57	4 -329 99 ± 67
	OTU or ASV assigned to	87	88	89	90
	genus and species (%)
	Reads assigned to plant	34	50	50	53
	pathogen (%)
	OTU or ASV assigned to	22	23	25	25
	plant pathogens (%)

Note: Results are presented for pooled OTUs/ASVs as indicated in the Materials and Methods.

Acknowledgements

We are very grateful to the RNSA analysts who provided the slides. We thank Julie Faivre d'Arcier for technical help on qPCR of E. alphitoides and Jean-Pierre Besancenot for comments on the manuscript. This research was financed by a grant from the French Ministry of Agriculture and Food through its Département de la santé des forêts (DSF) and benefitted from funding from the HOMED project funded under the Horizon 2020 research and innovation program under grant agreement No 771271. The ANSES, Laboratoire de la Santé des Végétaux, unité de mycology and the INRAE, UMR IAM are supported by a grant managed by the French National Research Agency (ANR) as part of the Investissements d'Avenir program (NR-11-L BX-0002-01, Laboratory of Excellence ARBRE)͘ The high throughput sequencing was conducted at the Genome Transcriptome Facility of Bordeaux (grants from ANR-10-EQPX-16). We are grateful to the INRAE MIGALE bioinformatics platform (https://migale.inra.fr/) for providing computational resources.