In the present study, the same water samples as in Katano et al. (2017) were used. Study sites are located in the headwater tributaries of the Ibo River in Hyogo, Japan (Fig. 1). Seven sites are located in the Akazai River (A1-A7) and 12 sites in the tributaries of the Ibo River (I1–I12). The orders of the streams and width of the tributaries are 1–2 and 3–5 m, respectively. A1–A7, I7 and I9, and I8 and I9 are located on the same tributary. At each sampling site, longitude, latitude, and altitude were recorded using a GPS device (Germin etrex 30×).

Figure 1.

Sampling sites and land use of the catchment area.

A 1 L sample was collected from the surface of the water at the central part of each stream in the study sites using a DNA-free polypropylene bottle. Samples were stored in a cooler box with a blank. The blank contained 1 L distilled water which was brought to the field and treated identically to the other water sampling bottles, except that it was not opened at the field sites. The sampling was conducted in late September, 2015. The collected water samples and blank were vacuum-filtered through 47 mm GF/F glass filters (pore size 0.7 μm, GE Healthcare, Little Chalfont, UK). The filter was then wrapped in commercial aluminum foil and stored at -20 °C before DNA extraction.

The protocols described by Uchii et al. (2016) were followed to extract DNA from the filters. In brief, each filter was incubated in 400 μl of Buffer AL (Qiagen, Hilden, Germany) and 40 μl of Proteinase K (Qiagen, Hilden, Germany). Then, 220 μl of TE buffer (10 mM Tris-HCl and 1 mM EDTA, pH: 8.0) was added to the filter and centrifuged. The dissolved DNA in the eluted solution was purified using a DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany) according to the manufacturer’s protocol.

For MiSeq sequencing, the fungal internal transcribed spacer 1 (ITS1) region of rDNA was amplified (Schoch et al. 2012). The methods of DNA analysis followed those described by Ushio et al. (2017) with some modifications, such as number of PCR cycles. The first-round PCR (first PCR) amplified the ITS1 region using the ITS1-F-KYO2 (5′- TAG AGG AAG TAA AAG TCG TAA -3′) and ITS2-KYO2 (5′- TAG AGG AAG TAA AAG TCG TAA -3′) primer set (Toju et al. 2012) fused with an Illumina sequencing primer and six random bases (N). The 1^st PCR was performed in a 12 μl volume with the buffer system of KODFX NEO (TOYOBO, Osaka, Japan), which contained 2.0 μl of template DNA, 0.2 μl of KOD FX NEO, 6.0 μl of 2× buffer, 2.4 μl of dNTP, and 0.7 μl each of the two primers (5 μM). The PCR conditions were as follows; an initial incubation for 2 min at 94 °C followed by 5 cycles of 10 s at 98 °C, 30 s at 68 °C for annealing and 30 s at 68 °C, 5 cycles of 10 s at 98 °C, 30 s at 65 °C and 30 s at 68 °C; 5 cycles of 10 s at 98 °C, 30 s at 62 °C and 30 s at 68 °C; 25 cycles of 10 s at 98 °C, 30 s at 59 °C and 30 s at 68 °C, and a final extension of 5 min at 68 °C. Eight replicate first-PCRs (per sample) were performed to mitigate the reaction-level PCR bias. Then, the duplicated first PCR amplicons (per sample) were combined, resulting in a template per sample for the second PCR. The PCR templates were purified using Agencourt AMPure XP (Beckman Coulter, Brea, California, USA) before the second PCR. The second PCR amplified the first PCR amplicons to add the MiSeqP5/P7 adapter and the 8X bases index sequences (Hamady et al. 2008) using the sequence primers. The second PCR was carried out with 12 μl reaction volume containing 1.0 μl of template, 6 μl of 2× KAPA HiFi HotStart ReadyMix (KAPA Biosystems, Wilmington, Washington, USA), 1.4 μl of each primer (2.5 μM), and 2.2 μl of sterilized distilled water. The PCR conditions were as follows; an initial incubation for 3 min at 95 °C followed by 12 cycles of 20 s at 98 °C, 15 s at 72 °C for annealing and extension, and a final extension of 5 min at 72 °C.

The second PCR amplicons were pooled and purified using AMPure XP. A target-sized DNA of the purified library (approximately 380–510 base pairs [bp]) was then excised using E-Gel SizeSelect (ThermoFisher Scientific, Waltham, MA, USA). The DNA sample was applied to the Illumina MiSeq platform at Ryukoku University, Japan. The sequence data were deposited in the Sequence Read Archive of the DNA Data Bank of Japan (accession number: DRA007786). See Appendix 1 for details of molecular analyses.

The procedures used for bioinformatics and data analyses followed those described previously (Matsuoka et al. 2016a; Ushio et al. 2017). The raw MiSeq data were converted into FASTQ files using the bcl2fastq program provided by Illumina. The FASTQ files were then demultiplexed using the commands implemented in Claident pipeline (Tanabe and Toju 2013; software available online: https://www.claident.org/). The forward and reverse sequences were then merged with each other. The total 611,092 reads (32,162 ± 10,784 reads per sample, mean ± SD, n = 19) were assembled using Claident v0.2.2016.07.05. First, short (<150 bp) reads were removed, then potentially chimeric sequences and sequencing errors were removed using UCHIME v4.2.40 (Edgar et al. 2011) and algorithms in CD-HIT-OTU (Li et al. 2012), respectively. The remaining sequences were assembled at a threshold similarity of 97% (Osono 2014), and the resulting consensus sequences represented molecular operational taxonomic units (OTUs). Then, singletons were removed. Through these procedures, 2,549 OTUs (602,304 reads) in total were obtained. An OTU table (i.e., matrix of OTUs and samples with sequence reads in each cell entry) was generated after this process. Then, coverage-based rarefaction was conducted to standardize the effect of sequencing depth (Chao and Jost 2012, rarefaction curves for each sample are shown in Suppl. material 1: Fig. S1). For each of the obtained OTUs, taxonomic identification was conducted based on the query-centric auto-k-nearest-neighbor (QCauto) method (Tanabe and Toju 2013) with NCBI database and subsequent taxonomic assignment with the lowest common ancestor algorithm (Huson et al. 2007) using Claident. The functional guild of each fungal OTU was estimated based on the FUNGuild database (Nguyen et al. 2016). Thirty-nine OTUs (718 reads) that were identified as non-fungal organisms were discarded and the remaining 1,956 OTUs (175,101 reads) were used in further analyses (Suppl. material 2: Table S1).

To test the effect of vegetation of the catchment area on fungal DNA assemblages, the gross area and vegetation of the catchment area of each sampling site was calculated by using geographic information system software (QGIS v.2.14.8). First, the catchment area of each site was calculated based on the GPS coordinates of the site. Then, by consulting the vegetation map made by the Ministry of the Environment of Japan, vegetation of the catchment area was divided into seven categories (deciduous broad-leaved forest, cedar-cypress plantation, ever-green broad-leaved forest, ever-green coniferous forest, grassland, cropland, and urban area) and the area of each category was calculated (Fig. 1).

The presence or absence of the OTUs was recorded for each sample regardless of the number of reads and used for all statistical analyses as binary data. All analyses were performed using R v.3.4.3 (R Core Team 2017). Note that in the present study, analyses described below were conducted on the datasets with unidentified OTUs (those that could not be identified even at the kingdom level). However, to investigate the effect of inclusion of unidentified OTUs, and the use of presence/absence data instead of read number, additional analyses were conducted by using datasets for only fungal OTUs (1,031 OTUs, see Suppl. material 2: Table S1), and by considering read counts of each OTU.

The variation of community composition of OTUs among sampling sites and the community dissimilarity of the OTUs among plots were visualized using non-metric multidimensional scaling (NMDS, Minchin 1987). Presence/absence data of OTUs for each site were converted into a dissimilarity matrix using the standardized effect size (SES) of the Jaccard dissimilarity index to show the community dissimilarity without dependence on richness gradient among sites, which is an issue of long-standing importance (Chase and Myers 2011). The SES was defined as: (β_obs - β_null) / β_sd, where β_obs is the observed β diversity, β_null is the mean of the null distribution of β diversity, and β_sd is the standard deviation of the null distribution. The null distribution was calculated based on 999 randomizations preserving both the site occurrence and the OTU richness. The randomized community data were generated using “randomizeMatrix” function in “picante” package. The SES value indicates the magnitude of the deviation from the expectation of a random assembly process, and positive and negative values indicate more and less β-diversity, respectively, than expected by chance. I12 is the only site located in an urban area, and the NMDS ordination indicates that the community in I12 is largely different from the others (Suppl. material 1: Fig. S2). Therefore, we excluded this site from further analyses. The correlation of NMDS structure with geographic coordinates (latitude and longitude) was tested through permutation tests (‘envfit’ command in the vegan package, 9999 permutations) to show the spatial structure of OTU composition among sites. Further, only the results for presence/absence data are shown. Then, to test whether the OTU compositions are more similar within a tributary than between tributaries, SES values were compared calculating 99% confidence interval (CI).

To test the effects of environmental variables (i.e., elevation, catchment area, and vegetation) and spatial distance on the OTU composition, a multivariate permutational analysis of variance (PERMANOVA) was implemented (‘adonis’ command in the vegan package, 9999 permutations). In these analyses, only one sample was used for each branch as samples from the same branch (A1-7 and I7-8) share most of their catchment. Therefore, eleven samples (A4, I1-6, and I8-11) were used for the following analyses and others were excluded. Preliminary analyses confirmed that the choice of sample inside the branch does not affect the results (data were used in Suppl. material 2: Table S1). Presence/absence data of ECM fungal OTUs for each plot were converted into a dissimilarity matrix using the modified Raup–Crick dissimilarity index (Chase et al. 2011). The Raup–Crick dissimilarity index is calculated based on the null model approach, which is akin to the SES value described above. Then, vegetation and spatial variables were calculated. Vegetation variables were extracted based on principal coordinate analysis (PCA). First, the vegetation data of the catchment area of each sampling site were transferred into the Euclidian distance matrix. Then, PCA vectors were extracted. The first two PCA axes explained the variations in the vegetation among sites of 86.9% and 11.2%, respectively, (Suppl. material 1: Fig. S3) and were used for the analysis. Spatial variables were extracted based on principal components of neighbor matrices (PCNM, Borcard et al. 2004). The PCNM analysis produced a set of orthogonal variables derived from the geographical coordinates of the locations of sampling sites. Two PCNM variables that best accounted for autocorrelation were used for the analysis.

After quality filtering and denoising, 602,304 reads of the original 611,092 reads were retained. After rarefaction analysis, 175,101 reads from 19 samples were grouped into 1,956 OTUs with 97% sequence similarity (Suppl. material 2: Table S1). In total, 770 OTUs were assigned as Ascomycota (38.4 % of the total number of fungal OTUs), 177 OTUs as Basidiomycota (8.8 %), and 38 OTUs as Chytridiomycota (1.9 %). The remaining 941 OTUs were not assigned to any phylum. The proportions of OTU numbers of each phylum were similar between the sampling sites except for I12 (Suppl. material 1: Fig. S4). At the order level, 552 Ascomycota OTUs belonged to 31 orders and 122 Basidiomycota OTUs belonged to 21 orders. The common ascomycete orders were Capnodiales (216 OTUs), Helotiales (88 OTUs), and Diaporthales (47 OTUs). The common basidiomycete orders were Tremellales (31 OTUs), Polyporales (15 OTUs), and Filobasidiales (12 OTUs). Thirty-six of 38 chytrid OTUs belonged to Rhizophydiales. At the genus level, 224 Ascomycota OTUs were assigned to 114 genera and 88 Basidiomycota OTUs to 55 genera. The OTU rich genera were Ramularia (36 OTUs) and Mycosphaerella (16 OTUs). The detailed taxonomic assignments are shown in Suppl. material 2: Table S1.

FUNGuild assigned 539 OTUs to functional guilds, 147 OTUs of which were saprotrophs and the others included pathogens (117 OTUs) and symbionts (19 OTUs) (Suppl. material 2: Table S1). The remaining 256 OTUs are those that were assigned to more than one functional guild, or whose function can change depending on circumstances. The detected OTUs included aquatic hyphomycetes, e.g., Flagellospora sp. (OTU_0413) and Tetracladium marchalianum (OTU_1491), that were dominantly detected from 19 and 16 out of 19 samples, respectively. On the other hand, OTUs belonging to terrestrial taxonomic or functional groups were also detected; plant decomposer (e.g., Mycena sp., OTU_0020); plant pathogen (e.g., Ciboria shiraiana, OTU_1523); and mycorrhizal fungi (e.g., Glomeromycetes sp. (OTU_1449) and Tuber sp. (OTU_0010)).

The Mantel test and NMDS ordination revealed biogeographic changes in the fungal OTU compositions. The dissimilarity of OTU composition among samples was significantly correlated with latitudinal and longitudinal distances (Mantel test, latitude, r = 0.585, P < 0.001; longitude, r = 0.510, P < 0.001). The NMDS ordination showed the geographical change of OTU composition among sites (Fig. 2, stress value = 0.154). The ordination was significantly correlated with the latitude and longitude (‘envfit’ function; latitude, r² = 0.737, P < 0.001; longitude, r² = 0.781, P < 0.001). The dissimilarity indices of OTU composition (SES value of Jaccard index) were significantly lower than 0 when compared between the sites in the same branch, higher when compared between distant branches, and intermediate when compared between adjacent branches (Fig. 3, based on the 99% confidence interval). This indicates that OTU compositions are significantly more similar inside the branch than expected by chance and more different between spatially distinct different branches. PERMANOVA showed that the differences in OTU composition among the branches were only related to the spatial variable (PCNM1) (F-model = 6.938, R² = 0.432, P = 0.009) and did not significantly correlate with the other environmental variables including elevation and vegetation PCA axis (Table 1). In addition, using datasets composed only by fungal OTUs and taking into account reads abundance of each out, it did not affect the results. The results for these additional analyses are provided in Appendix 2.

Figure 2.

Dissimilarity of the DNA assemblages among sites as revealed by nonmetric multidimensional scaling (NMDS) ordination (stress value = 0.154). Numbers are identical with site numbers in Fig. 1. Fungal DNA assemblages were similar within the same branch (A1-7) and changed with latitude and longitude across the study site.

Figure 3.

Standard effect size (SES) of the Jaccard index for comparison of the compositional dissimilarity within a same branch and between the adjacent and distant branches. Error bars indicate 99% confidence interval. The horizontal dotted line represents SES = 0, indicating a non-significant effect.

River water contained both phylogenetically and functionally diverse fungal DNA. In the present study, phylogenetically diverse fungi, primarily ascomycetes and basidiomycetes, were detected in the river, as shown in the recent studies using molecular methods (Duarte et al. 2015, Wurzbacher et al. 2016; LeBrun et al. 2018). In river water, aquatic fungi have diverse lineages and functional groups, including saprophytes, which decompose plant substrates, and parasites of plants and animals. In the present study, some genera known as aquatic hyphomycetes (e.g., Flagellospora and Tetracladium; Raja et al. 2018; Grossart et al. 2019) were detected. The OTUs belonging to these taxa are thought to degrade plant substrates that have entered the river from land (fallen leaves and wood). In particular, OTUs such as OTU_0413 (Flagellospora sp.) and OTU_1491 (Tetracladium marchalianum) were detected in almost all samples, suggesting that they are widely distributed in the surveyed landscape. These aquatic fungi can be detected based on the eDNA obtained from the river water owing to the presence of spores or mycelial fragments in the water.

OTUs of terrestrial fungi such as the plant decomposer Mycena and the mycorrhizal fungi Tuber and Glomeromycota were detected, as well as the OTUs of aquatic fungi. DNA of these fungi may be released from the propagules, which flowed into the rivers from land. For example, Basidiomycota acting as plant decomposers often produce epigeous sporocarps and disperse their spores into air. Such airborne spores have the potential to enter the river system. On the other hand, species in the genus Tuber produce hypogenous sporocarps and do not disperse their spores into air. This indicates that the propagules and mycelial fragments in soil and plant materials (e.g., decomposing leaves) flow into the river via the surface current or with the direct addition of these substrates (Voronin 2014). Recent studies that investigated the environmental fungal DNA of lakes and rivers also detected a significant amount of DNA for terrestrial fungi, as well as for aquatic fungi (Deiner et al. 2016, Khomich et al. 2017; LeBrun et al. 2018). Our results are consistent with the results of these previous studies.

The fungal OTU compositions in the studied landscape were more similar within the same branch than between branches. This may be because some fungal DNA, which entered upstream flowed downstream and was therefore detected both upstream and downstream in the branches. Indeed, recent evidence suggests that rivers work as a dispersal pathway of terrestrial fungi (Deiner et al. 2016; LeBrun et al. 2018). Other possibilities are that the same aquatic fungal species live across the branches, or that the same terrestrial fungal DNA enters the river at several different sites. For example, Tsui et al., (2001) investigated aquatic fungi in upstream and downstream of the same river and reported the presence of common species in both areas of the stream. However, there is currently no information in support of either of these possibilities. Further investigations into the composition and distribution of aquatic fungi along the river, the composition of fungal DNA flowing downstream, as well as the air and surface currents may help to clarify the mechanisms behind these results.

Fungal DNA assemblages show geographical structures at a scale of around 15 km in the forest site. This means that fungal DNA assemblages are more similar between spatially closer branches than between distant ones. Such geographical structures have also been observed by Khomich et al. (2017); fungal DNA assemblages along latitudes and longitudes were investigated in 77 lakes in the Scandinavia Peninsula. In the lakes, the surrounding rivers brought and accumulated the DNA of terrestrial fungi from the surrounding areas (Khomich et al. 2017). Our study shows that such geographical structures can also be observed in streams, where water has a directional flow and is constantly replaced. The geographical structure of fungal DNA may be generated by both aquatic and terrestrial fungal DNA. For aquatic fungi, environmental differences (e.g. water chemical properties or current velocity/discharge) across branches or dispersal limitations can affect the fungal communities (Heino et al. 2014; Duarte et al. 2017), causing geographically closer branches to share more fungal species. Similarly, environmental factors such as vegetation type, soil chemical properties, and dispersal limitation contribute to the formation of similar terrestrial fungal communities in the same vicinity (Toljander et al. 2006; Kadowaki et al. 2014). The propagules of these fungi flow into the nearby river, causing the DNA assemblages of geographically closer branches to resemble one another.

In the present study, geographic variations in the fungal DNA assemblage were explained only by spatial variables, and no significant relationship between fungal DNA assemblage and the vegetation of the catchment area was detected. These results suggest that in the study landscape, the spatial structure of the fungal DNA assemblages may be influenced by spatial factors such as dispersion limitation (Talbot et al. 2014; Peay et al. 2016) or unmeasured environmental factors. Previous studies have shown that the community composition of aquatic fungi is more strongly influenced by spatial variations in the environment, such as water quality (e.g., pH), compared to spatial distance among sites (Heino et al. 2014; Duarte et al. 2017). For example, LeBrun et al. (2018) investigated fungal DNA assemblage in river water with P enrichment problems and showed that the fungal DNA assemblage changes with P concentration. In our study sites, no such nutrient enrichment problems have been reported to date (Ministry of Land, Infrastructure, Transport and Tourism 2007) and the range of variation in the environment is considered to be smaller than that of previous studies (Heino et al. 2014; Duarte et al. 2017), which dealt with community changes at regional and continental scales. However, it is highly likely that unmeasured water qualities influence DNA assemblages in our study sites. In addition, the discharge of the surveyed river may affect the DNA compositions via changing the size of particles carried in the water column. In the present study, the brown particles, such as the fine particle organic matter (FPOM), were attached to the filter and fungal DNA attached to the FPOM might have been detected. Wurzbacher et al. (2016) showed that the detected fungal DNA can vary depending on the size and amount of FPOM. These unmeasured factors and relative importance of environmental variables and spatial distance should be considered in future studies. In a forest ecosystem, the types of vegetation present have a stronger effect on fungal communities compared to other environmental variables (Peay et al. 2013). Therefore, if the proportion of terrestrial fungal DNA is high, the fungal DNA assemblage may be explained by the vegetation present. The lack of a relationship between the vegetative and fungal DNA assemblages can be partly attributable to the fact that the catchment area and distribution of fungal individuals of potential DNA origin did not correspond. For example, at a landscape scale with tens of kilometers, similar to our study site, fungal DNA would possibly flow or be blown from surrounding forests, including outside the catchment area. Otherwise, fungal DNA can be derived from areas closer to the river catchment area. There is limited information available on the dispersal ability of fungi. However, a previous study has reported that for terrestrial fungi, most spores fall within several meters of the sporocarps (Li et al. 2005; Galante et al. 2011). This indicates that spores of fungal species and individuals colonizing more than several meters apart from the river should have little chance to fall into the river. Further studies are required to specify the range from which most fungal DNA is derived.

In the present study, the ITS region was used as a genetic marker, but many OTUs remained unidentified. Compared to other terrestrial habitats, fewer genomic studies have been conducted in aquatic habitats. As a result, fewer ITS sequence records are available and new phylogenetic groups are expected to be found (Grossart et al. 2019; Khomich et al. 2018). Furthermore, for basal fungi living in an aquatic habitat, such as chytrids, more conservative regions such as large subunit (LSU) and small subunit (SSU) regions of rDNA are often used (Nilsson et al. 2019). However, few studies have investigated whether a difference in genetic markers affects the analytical results (Nilsson et al. 2019). Therefore, this may be an interesting topic for future studies to investigate. This database deficiency and handling of resulting unidentified OTUs as well as the selection of appropriate marker genes and primer limitation/specificity towards a particular group can be a key problem in fungal DNA metabarcoding.