Five wetlands from Central Italy were considered in this study: (1) Colfiorito Marsh, in Umbria Region (43°01′23″ N, 12°52′36″ E) at 756 m a.s.l.; (2) Lake Trasimeno, in Umbria Region (43°08′05.5″ N, 12°06′04.6″ E) at an average altitude of 257 m a.s.l.; (3) Lake Vico in Lazio Region (42°18′58.40″ N, 12°10′5.89″ E), 507 m a.s.l.; (4) Lake Chiusi in Tuscany Region (43°03′22.11″ N, 11°57′55.79″ E) at 252 m a.s.l.; (5) the wetland of “Le Morette” (43°48′30.38″ N, 10°48′20.14″ E) at Padule of Fucecchio Marsh, in Tuscany, at 13 m a.s.l. (Figure 1). These five ecosystems were chosen to cover a large variability in terms of: (i) surface area, from large lakes such as Trasimeno and Vico, to the small wetland of Le Morette and Colfiorito marshes; (ii) water depth, from the deep waters of Vico, to the shallow waters of Le Morette or Colfiorito marshes; (iii) microclimate, from temperate Apennine areas such as Colfiorito, to the submediterranean climate of Vico or Chiusi lakes; (iv) previous reed retreat documented by aerial images (Gigante and Venanzoni, 2012; Gigante et al., 2014).

The sampling design was based on the results of former investigations in the same areas (Lastrucci et al., 2016b). Two types of stands could be identified for sampling: healthy stands (group N on Supplementary Table S1) and RDBS stands (group Y on Supplementary Table S1), the latter showing a number of RDBS diagnostic traits and, above all, the clumping habit. This is a clustered growth of the reed stems deriving from an uncontrolled outgrowth of dormant buds caused by the breaking of apical dominance, and represents a widely acknowledged and easily detectable trait of reed die-back (Armstrong et al., 1996; van der Putten, 1997; Clevering, 1998; Dinka and Szeglet, 2001; Gigante et al., 2011). Sampling was carried out during two different periods: in February 2015, bulk and rhizosphere soils were sampled from three healthy stands and three RDBS stands from four wetlands (Colfiorito, Trasimeno, Vico, and Chiusi) (48 samples); in September 2015, two healthy stands and two RDBS stands were sampled from all five wetlands, collecting, where possible, bulk, rhizosphere, and roots (48 samples). Bulk soil (40 g for molecular analysis and 300 g for chemical analysis) was collected at an approximate depth of 30 cm, 50 cm away from the plant, using a soil core sampler. Belowground part of the stem with roots were collected; roots were cut and thoroughly washed in distilled water for several times, while the soil attached to the roots was considered as rhizosphere. Soil and root samples were stored at -20°C until DNA extraction. For a summary of the samples used in this study, see Supplementary Table S1.

Two different procedures were used for DNA extraction, one for soil samples (bulk soil and rhizosphere) and the second one for root samples. Soil samples weighing 40 g were mixed, freeze dried for 48 h, and ground; from these samples, 250 mg were used for total soil DNA extraction using the PowerLyzer^TM PowerSoil^® DNA Isolation Kit (Mo Bio Laboratories, Carlsbad, CA, United States) according to the manufacturer’s instructions, except that samples were homogenized in a Geno/Grinder 2000 (SPEX CertiPrep, Metuchen, NJ, United States) at 1500 rpm for three times, 30 s each. For root samples, 100 mg of fresh material were ground in a mortar with liquid nitrogen; DNA was extracted using the DNeasy Plant Mini kit (Quiagen GmbH, Hilden, Germany) according to the manufacturer’s instructions.

The internal transcribed spacer (ITS) 1 region of ribosomal DNA was used as a marker. Each sample was barcoded with 10-nucleotide multiplex identifier (MID) primer tags before pooling PCR products (Supplementary Table S1). Amplification strategy was as described by Sapkota and Nicolaisen (2015): the PCR reaction mixture consisted of 1× PCR reaction buffer, 1.5 mM of MgCl₂, 0.2 mM of dNTPs, 1 μM of each primer, and 1 U of GoTaq Flexi polymerase (Promega Corporation, Madison, WI, United States) in a total volume of 25 μl containing 24 μl of reaction mixture and 1 μl of template. Tagged PCR amplicons from the 96 samples were pooled in equimolar amounts, electrophoresed on a 1.5% agarose gel and the visible smear of 200–400 base pairs was excised and purified from the gel using QIAquick Gel Extraction Kit (Quiagen, Hilden, Germany). The pooled samples were sequenced by Eurofins MWG (Ebersberg, Germany) on a GS Junior 454 Sequencer (Roche Diagnostics) using titanium chemistry. The raw data were deposited on the Sequence Read Archive website¹ under the SRA number SRP111776.

Sequences were analyzed using the bioinformatics pipeline available in QIIME v. 1.8 (Caporaso et al., 2010). De-multiplexing and quality filtering was carried out using split_libraries.py command using quality score of 50, discarding reads less than 150, homopolymer length of 10 nucleotides, maximum number of ambiguous bases of 6 and other default settings. The ITS1 region was extracted using ITSx extractor version 1.0.6 (Bengtsson-Palme et al., 2013). ITS1 reads were then clustered using the pick _otus.py script at 98% similarity level (Vettraino et al., 2012) using usearch61 with a minimum cluster size of two, thus excluding singletons (Edgar, 2010). For taxonomic assignment, a reference database described earlier (Sapkota and Nicolaisen, 2015) was used. In addition, for species identification, a representative sequence from each of the most abundant operational taxonomic unit (OTUs) with at least 100 reads in total was further subjected to Basic Local Alignment Search Tool (BLAST) matches at NCBI. At least the top five BLAST hits were evaluated: OTUs with identity and coverage ≥ 95% were assigned at species level, or at genus or family level if lower. Only reads assigned to oomycetes were retained in the OTU table.

The OTU table was used for diversity based calculations. Data visualization and statistical analysis were carried out in R (R Core Team, 2015). A species accumulation curve was generated using the specaccum function available in the ‘vegan’ package (Oksanen et al., 2015). Samples representing less than 400 reads were removed before performing diversity based calculations. Species richness and evenness were calculated by rarifying the OTU table at an even sampling depth of 400 reads, while the OTU table was transformed to relative abundance before beta diversity based calculations. The species evenness was evaluated by Pielou’s index and the species richness was estimated as the number of OTUs. β-diversity distance matrices were obtained using Bray–Curtis dissimilarity and visualized using principal coordinates analysis (PCoA) plot. Permutational multivariate analysis of variance (Anderson, 2001) with a permutation number of 999 was carried out to test the null hypothesis of no differences among a priori defined groups using the adonis function in ‘vegan’ package.

Covariance biplots were calculated using composition data analysis on centered log-ratio transformed OTU table (Gloor and Reid, 2016), removing samples less than 400 reads.

The OTU table was also used to perform the Dufrene-Legendre indicator species analysis in R, using the labdsv package (Roberts, 2012) to test the OTUs associated with die-back. OTUs with significant values (P < 0.05) and an indicator value > 0.25 were set as cut-off values to define indicator species (Dufrene and Legendre, 1997). The difference in relative abundance of OTUs was determined by non-parametric statistical method using Kruskal–Wallis test for ANOVA and Wilcoxon rank sum tests for pairwise comparison. The network analysis was made in QIIME using OTU table (with the command make_otu_network.py) and the resulting network files were visualized in cytoscape². Samples and OTUs were subjected to network analysis via the spring-embedded algorithm using eweights, where samples represent larger nodes than OTUs. All OTUs are connected to the sample to which they belong with a line (edge). Shared OTUs are thus connected to several samples, and occupy the center of network, whereas OTUs that only occur in one or few samples occupy the periphery of the network. Samples are thus clustered based on shared OTU’s and the degree to which samples cluster is based on the number of OTUs shared between samples.

Using the indicator species analysis, 12 OTUs were found to have significantly higher abundance in die-back samples in the total dataset; they were aligned together with the ITS1 sequences of the accessions of Pythiogeton available in GenBank. ITS1 sequences of the two most abundant Pythium species in our dataset were used as outgroup. The Neighbor-joining tree was assembled with MEGA 7.0 (Kumar et al., 2015) and tested with a bootstrap of 1000 replicates, to ascertain the reliability of any given branch pattern.

Concentrations of heavy metals (Cd, Pb, Zn, Cr, Ni, Cu, Al) in soil samples and in plant tissues were determined by Inductively Coupled Plasma Optical Emission Spectrometry (ICP-OES, Ultima 2, HORIBA Scientific) equipped with an ultrasonic nebulizer (U-5000AT, CETAC Technologies). The nebulizer gas flowrate was kept at 1 L min^-1 and the plasma operated at a radiofrequency power of 1000 W. Supplementary Table S2 shows the wavelength used. Instrumental detection limits were in the range of 0.1–1.5 mg kg^-1. Commercially produced (ICP multi-element standard solution IV CertiPUR^®, VWR Merck Chemicals and Reagents) standard solutions (1000 mg L^-1) in nitric acid were used to prepare appropriate elemental calibration standards. An acid digestion was performed on the soil (Ministero dell’Ambiente e della Tutela del Territorio [MATT] and Agenzia per la Protezione dell’Ambiente e per i Servizi Tecnici [APAT], 2005) and on the plant tissue samples (Du Laing et al., 2003) prior the analytical determination. The soils were air-dried, disaggregated using a mortar and pestle to pass through a 2 mm mesh sieve, dried at 105°C for 24 h and digested as follows: 8 ml of ultrapure nitric acid (Millipore Suprapur^®, 65%) and 2 ml of ultrapure solution of hydrogen peroxide (Millipore Suprapur^®, 30%) were added to 0.200 g of soil samples and digested in a Mars Microwave Oven, working at a power of 1000 W. Microwave digestion consisted of two steps: 130°C (200 psi) for 1 min, 180°C (300 psi) for 10 min. The mixture was cooled, filtered (Whatman Grade No. 42, particle retention 2.5 mm) and diluted with ultrapure water to 50 ml. Plant samples were preliminarily dissected into roots and rhizomes and dried at 70°C for 48 h. Acid digestion of tissues (0.4 g) was carried by the same procedure used for soils, but ultrapure water was added to reach the volume of 25 ml. The measure of the heavy metal was based on two replicates of the overall analytical procedure for each sample of soil and plant tissue examined.

In total, 235065 reads were obtained from 454-pyrosequencing. After quality control, ITS extraction and taxonomic assignments, a total of 207661 reads were clustered into 523 OTUs assigned to oomycetes, excluding singletons. Two samples failed to yield any reads; thus, reads were obtained from 94 samples in total. The number of reads per sample revealed high variation (min, mean, max: 4, 2464, 13461).

The accumulation curve showed that the number of OTUs hardly approached a plateau phase (Figure 2), indicating a higher oomycetes diversity in the area than was covered in this study and that diversity would increase with increasing sample number. In our dataset, Pythiales dominated the reads (66%), followed by Lagenidiales (25%), Saprolegniales (6.5%) and Peronosporales (2.4%), while at genus level we detected the presence of three dominating taxa: Pythiogeton (30.9%), Pythium (29.9%), and Lagena (22%) (Figure 3).

The incidence (the number of samples in which a certain OTU occurred) was checked for the 523 OTUs identified in this study: 166 OTUs were found in only one sample, while only one OTU (denovo0; Lagena) was found in more than 50% of the samples (see Supplementary Table S3). Lagena sp. (denovo0) had a significantly higher relative abundance in root samples in comparison to rhizosphere (Kruskal–Wallis test, P = 0.023) and bulk soil samples (Kruskal–Wallis test, P = 0.017) (Figure 4). The species richness was compared in relation to different sites, plant compartments (root, rhizosphere, or bulk soil), temporary or permanent flooding (i.e., absence or presence of die-back). Interestingly, none of the studied factors significantly influenced OTUs richness. The species evenness estimated by Pielou’s index, is significantly influenced by the compartment (df = 2; P = 0.017). Pairwise Wilcoxon rank sum tests revealed that the evenness is significantly different between root and bulk soil samples (P = 0.018).

The effects of the geographical location, the compartment (bulk soil, rhizosphere soil, or root sample), die-back occurrence, and the available chemical data were subjected to permutation multivariate analysis of variance to evaluate factors that affected the structure (beta diversity) of oomycetes communities. The factors that were found to be significant in explaining the diversity were die-back occurrence (R² = 0.049, P = 0.001), compartment (R² = 0.043, P = 0.009), and the interaction between lake and die-back occurrence (R² = 0.073, P = 0.011) (Table 1). Among the heavy metals, only cadmium concentration (mg⋅Kg^-1) is significantly related to oomycetes distribution (R² = 0.035, P = 0.002) (Table 2).

In support of this, a PCoA plot revealed clustering based on the plant health status, demonstrating the importance of die-back in explaining the diversity of oomycetes (Figure 5). The clustering of samples based on die-back was supported by covariance biplots (Supplementary Figure S1). The influence of die-back in oomycetes communities was also visible via network based analysis (Figure 6). The samples with higher number of shared OTUs were closer and clustered in the center of the network. The majority of die-back samples clustered in the lower core of the network. The network analysis also showed that a large part of OTUs were shared between healthy and RDBS samples.

Indicator species analysis of the total dataset revealed 12 OTUs that were associated with RDBS (Table 3).

Those 12 OTUs belonged to the genus Pythiogeton. The relative abundance of Pythiogeton was significantly higher in die-back samples in comparison to healthy (P = 0.0002) (Supplementary Figure S2). A phylogenetic tree (Figure 7) shows that eight of the OTUs (denovo594, denovo1113, denovo208, denovo8, denovo3, denovo245, denovo302, denovo246) are similar to Pythiogeton microzoosporum, although the clades do not show a high bootstrap support. Denovo1 and denovo284 are similar (bootstrap = 98) and close to the clade containing Pythiogeton oblongilobum (bootstrap = 63). Denovo34 and denovo152 cluster together (bootstrap = 100) and they are similar to Pythiogeton abundance.