Temporal variation in community structure of zoosporic fungi in Lake Biwa, Japan

Zoosporic fungi play an important role in aquatic environments, but their diversity, especially that of parasitic fungi of phytoplankton, has still not been fully revealed. We conducted monthly analyses of the community structure of zoosporic fungi at a pelagic site in Lake Biwa, Japan, from May to December 2016. Metabarcoding analysis, targeted to a large subunit region of ribosomal DNA in the nano-size fraction of particles (2−20 μm), was carried out on the samples. We also counted large phytoplankton and chytrid sporangia attached to the hosts. We detected 3 zoosporic fungal phyla (Blastocladiomycota, Chytridiomycota and Cryptomycota) within 107 operational taxonomic units (OTUs), in which Chytridiomycota was the most diverse and abundant phylum. Few fungal OTUs overlapped between months, and specific communities were detected in each month. These results showed that diverse zoosporic fungi with high temporal variability inhabited the lake. Five large phytoplankton species were found to be infected by chytrids: Staurastrum dorsidentiferum, S. rotula, Closterium aciculare, Asterionella formosa and Aulacoseira granulata. Some chytrids were detected by metabarcoding analysis: Zygophlyctis asterionellae infecting A. formosa, Staurastromyces oculus infecting S. dorsidentiferum and Pendulichytrium sphaericum infecting A. granulata. One OTU detected in association with infected C. aciculare by microscopic counting might have been an obligate parasitic chytrid of C. aciculare. The results indicated that a combination of metabarcoding and microscopic analysis revealed more information on zoosporic fungi, including those that are parasitic.

Recently, the taxonomy of zoosporic fungi has been greatly expanded with the development of molecular phylogenetics and ultrastructure analysis (Frenken et al. 2017). However, owing to missing DNA reference data, parasitic fungi are often not identifiable, even by deep amplicon sequencing. Phylogenetic analysis could assign taxonomic positions of unidentifiable sequences and shed light on their ecology based on phylogenetic positions (Seto et al. 2020). Counting fungal sporangia attached to host phytoplankton cells is helpful in identifying the parasitic fungi of phytoplankton. Combining microscopic counts and metabarcoding analysis may improve the detection of parasitic fungi and allow us to understand the community structure of zoosporic fungi.
In Lake Biwa, the largest lake in Japan, parasitic chytrids have been observed to infect the dominant green alga, Staurastrum dorsidentiferum, and to influence its temporal dynamics and material flows (Kagami & Urabe 2002. Diverse fungi including Cryptomycota and Aphelidiomycota, as well as Chytridiomycota, have also been detected from different phytoplankton species in Lake Biwa by analysis of the single sporangium-like cell attached to a phytoplankton cell or colony (Ishida et al. 2015). The spatio-temporal dynamics of these host− parasite relationships in Lake Biwa remain uncertain. High-throughput sequencing of the internal transcribed spacer 2 (ITS2) region of ribosomal DNA (rDNA) revealed fungal communities in Lake Biwa; however, only 1 zoosporic fungal phylum (Chy tridiomycota) was detected (Song et al. 2018). The ITS2 region of chytrids is too variable for alignment by phylogenetic analysis, and not many chytrid ITS2 sequences have been deposited in databases yet (Frenken et al. 2017, Heeger et al. 2019. The large subunit (LSU) region of rDNA has less genetic variation than the ITS2 region, but has adequate variation for chytrids, and many sequences of known chytrid species have already been deposited in databases and have been applied to analyze chytrid communities in several lakes (Lefèvre et al. 2012, Wurzbacher et al. 2016). In the present study, we determined the community structure of zoosporic fungi using highthroughput sequencing and phylogenetic analysis of the LSU region of rDNA extracted from nano-sized particles (2−20 μm), to clarify the diversity and tem-poral variation in the zoosporic fungal community in Lake Biwa. Chytrid sporangia attached to phytoplankton were also counted by staining, using epifluorescence microscopy to identify the parasitic chytrids infecting host phytoplankton.

Sample collection
Water samples were collected monthly from 6 May to 12 December 2016, when phytoplankton were abundant, with a Van Dorn sampler (volume: 10 l) from each depth at 1, 3, 5, 10 and 15 m at Stn 3 (70 m deep, 35° 19.8' N, 136° 08.9' E) in the north basin of Lake Biwa (mean depth, 41.2 m), Japan. For molecular analyses, equal volumes (500 ml) of water collected from each depth were mixed and combined into 1 sample; this was because chytrid distribution and any environmental parameters were vertically homogeneous above 15 m during the study period (Figs. S1 & S2 in the Supplement at www.int-res. com/articles/suppl/a087p017_supp.pdf). All biotic and abiotic parameters were therefore calculated as averages throughout the 0−15 m water column (Table S1). The samples for molecular analysis were transferred to the laboratory within 2 h in an insulated container to keep the samples cold. Aliquots of the samples from each depth for counting phytoplankton and chytrid sporangia cells attached to the host phytoplankton were preserved immediately after collection with buffered formalin at a final concentration of 1%. Vertical profiles of water temperature were determined with a CTD profiler (JFE Advantech, ACL1183-PDK) on each sampling occasion.

Phytoplankton and sporangia counts
Large phytoplankton (> 20 μm) species were identified and counted under a light microscope (Olympus, BX-50) at a magnification of 400×. Chytrid sporangia attached to the host cells were stained with calcofluor white and counted with blue fluorescence under excitation and emission wavelengths of 385 and 440 nm, respectively (Kagami et al. 2007a), using an epi-fluorescence microscope (Olympus, BX-51) at a magnification of 100−1000×. The prevalence of infection for each phytoplankton species was calculated as the percentage of phytoplankton cells with chytrid sporangia out of the total number of phytoplankton cells.

DNA extraction
The combined water sample from each depth was pre-filtered through a 20 μm mesh net for separating zoosporic fungi (2−5 μm in diameter) from other, larger fungi and those attached to the host phytoplankton and substrates. Next, 500 ml of the filtrate were filtered again through a 2.0 μm pore size membrane filter (Millipore) to collect the zoospores. DNA was extracted directly from the membrane filter with the DNeasy ® Plant Mini Kit (Qiagen ® ). DNA concentration was measured with a Qubit fluorometer (Invitrogen). The extracted DNA was frozen at −20°C until analysis.

Metabarcoding analysis using Illumina MiSeq
Metabarcoding analysis targeted the fungal-specific LSU regions, which were amplified using the primers LR3R (5'-GTC TTG AAA CAC GGA CC-3') (Rehner & Samuels 1995) and LR5 (5'-TCC TGA GGG AAA CTT CG-3') (Hopple & Vilgalys 1994). All PCR reactions were conducted in 20 μl mixes with Takara Ex TaqHS mix, 0.2 μM of both forward and reverse primers and ~1 ng of extracted DNA. The first PCR cycling consisted of initial denaturation at 94°C for 2 min; followed by 25 cycles of denaturation at 94°C for 30 s, annealing at 52°C for 30 s and elongation at 72°C for 30 s; and then a final extension at 72°C for 5 min. To attach the index sequence to the first PCR products, the second PCR products were processed with the primer set 2 nd F: AAT GAT ACG GCG ACC ACC GAG ATC TAC AC-Index2 (8 bp)-ACA CTC TTT CCC TAC ACG ACG C and 2 nd R: CAA GCA GAA GAC GGC ATA CGA GAT-Index1 (8 bp)-GTG ACT GGA GTT CAG ACG TGT G (Table S2). The second PCR cycling consisted of initial denaturation at 94°C for 2 min; followed by 8 cycles of denaturation at 94°C for 30 s, annealing at 60°C for 30 s and elongation at 72°C for 30 s; and then a final extension at 72°C for 5 min. The library quality was assessed using a Fragment Analyzer system. The libraries were sequenced on an Illumina MiSeq platform and 250 bp paired-end reads were generated. We obtained 1 112 173 raw reads from all samples collected, and all raw reads were deposited in the NCBI Sequence Read Archive (SRA) database under accession number SRP119396.

DNA amplification and sequencing
Barcodes and PCR primers were removed from the raw reads by the fastq_barcode_splitter from the Fastx toolkit software (ver. 0.0.13.2). Low-quality reads (values < 20) were trimmed using Sickle software (ver. 1.33). High-quality reads were processed to merge paired-end reads of each sample through FLASH software (ver. 1.2.10) (Magoč & Salzberg 2011). After these processes, the reads obtained were screened by USEARCH software (ver. 8.0.1623-i86linux64) to find chimera sequences (Edgar 2010). Finally, effective reads were generated by removing these chimera sequences, and singletons were removed. UPARSE software (Edgar 2013) was used to cluster effective reads of all samples. Effective reads with 97% similarity were classified into 1 operational taxonomic unit (OTU). Taxonomic assignments were based on assign_taxonomy.py scripts of Quantitative Insights Into Microbial Ecology software (QIIME, ver. 1.9.0) (Caporaso et al. 2010). Data that contained the term 'no BLAST hit' in the taxonomic annotation were removed. All OTUs that were not assigned to fungi were discarded from the analysis. OTU sequences belonging to fungi were submitted to the DNA data bank of Japan (DDBJ) Mass Submission System (MSS) under accession numbers LC323648−LC324677. Zoosporic fungal sequences (i.e. phyla Blastocladiomycota, Chytridiomycota, Cryptomycota and Neocallimastigomycota) were selected for the downstream analysis.

Statistical analysis
OTU composition and relative abundances of zoosporic fungal assemblages derived from the DNA amplification and sequencing results in each month were calculated to clarify taxonomic diversity in the zoosporic fungal community in the lake. In order to visualize the number of unique OTUs detected in each month, we used Gephi software (ver. 0.9.1), which is a tool for discovering the connectivity among the OTUs tested. Prior to the statistical analyses, a rarefied OTU table was generated by QIIME rarefaction analysis script with the rarefied depth in the minimum sequence number. Chao-1 (Chao 1984) and Shannon's H (Shannon 2001) indices were calculated as representatives of alpha diversity in fungal communities.

Phylogenetic analysis
In order to determine which OTUs were closely related to those of cultured chytrids by molecular phylogenetic analysis, we created a data set of LSU rDNA sequences of the class Chytridiomycetes. Abundant and/or unique Chytridiomycetes OTUs were selected and added into the data set. Three species of the class Monoblepharidomycetes were selected as outgroup taxa. Sequences were automatically aligned with MAFFT (Katoh & Standley 2013). Ambiguously aligned regions were excluded using 'trimAl' (Capella-Gutiérrez et al. 2009) with the 'gappyout' option. The maximum likelihood tree was inferred using RAxML 8.2.7 (Stamatakis 2014). We ran an analysis under the GTRGAMMAI model and used the '-fa' option to conduct a rapid bootstrap analysis with 1000 replicates combining 200 searches for the optimal tree.

Physicochemical environment and host phytoplankton
According to vertical profiles of water temperature, lake water was vertically stratified from June to September, when the temperature gradually decreased below 5 m (Fig. S1). However, most of the physicochemical parameters other than temperature were vertically homogeneous above 15 m even in the stagnation period (Fig. S1).
We found a total of 12 large phytoplankton species comprising 1 dinoflagellate (Ceratium hirundi nella), 6 green algae (Closterium aciculare, Micrasterias hardyi, Staurastrum arctiscon, S. dorsidentiferum, S. rotula, Xanthidium hastiferum) and 5 diatom species (Aulacoseira granulata, Asterionella formosa, Aulaco seira nipponica, Fragilaria crotonensis, Synedra sp.) (Fig. 1). According to the vertical profiles of the dominant phytoplankton species and chytrid sporangia that attached to their hosts, all were almost uniformly distributed throughout the 0−15 m water column (Fig. S2). Therefore, average abundances throughout the 0−15 m water column are shown in the results. Green algae were always the dominant group at > 50% of the total counted phytoplankton abundance, while diatoms and dinoflagellates comprised < 38 and <12%, respectively, throughout the study period. S. dorsidentiferum was the most dominant species among the phytoplankton (> 30% of total abundance) from May until September and had its highest density of 18.0 × 10 3 cells l −1 in August, while M. hardyi increased abruptly from October and had its highest density of 78.0 × 10 3 cells l −1 in December.
Chytrid sporangia were found on 5 out of 12 large phytoplankton species (Figs. 2 & 3). The prevalence of infection differed among phytoplankton species: the percentages of infected S. dorsidentiferum 20 Fig. 1. Species composition and cell densities of phytoplankton from 6 May to 12 December 2016 at Stn 3 in the north basin of Lake Biwa, Japan ranged from 4.7 to 25.2%, while chytrid sporangia appeared throughout the sampling periods (Fig. 2). The prevalence of infection in S. rotula ranged from 50 to 100% from September to November. Chytrids on C. aciculare appeared only in June and July, when 3.7 and 23.9% of the total abundance, respectively, were infected. Chytrids were found on A. formosa in May, October and November, and the prevalence of infection was relatively high, between 50 and 95.2%. A. granulata was mainly infected in October and November, at 5.2−24.6%. The total density of chytrid sporangia was highest in October with 7.6 × 10 3 cells ml −1 , and ranged between 0.3 × 10 3 and 2.4 × 10 3 cells ml −1 in other months (Fig. 2). No infected cells of M. hardyi were found, even when it dominated in November and December.

Sequence data
After filtering, 956 761 reads included a total of 13 757 OTUs in the final matrix. A total of 22 264 reads and 311 OTUs were assigned as zoosporic fungal sequences. An additional 204 OTUs were removed from the 311 OTUs as singletons, leaving 107 zoosporic fungal OTUs as the final data set, comprising 22 060 reads. OTUs assigned to the phylum Neocallimastigomycota were excluded because phylogenetic analysis indicated that these sequences were not related to Neocallimastigomycota (see Section 4); as a result, 101 OTUs and 15 516 reads were used for the analyses (Table 1).

Zoosporic fungal diversity and taxonomic characteristics in Lake Biwa
Three zoosporic fungal phyla, namely Blastocladiomycota (1.5% of reads and 2.0% of OTUs), Chytridiomycota (97.4 and 93.0%) and Cryptomycota (1.1 and 5.0%), were detected. Chytridiomycota was the most diverse and abundant phylum among them. The dominant orders were Chytridiales (38.6% of reads and 27.7% of OTUs), Rhizophydiales (32.5 and 28.7%) and Spizellomycetales (10.4 and 14.9%). The A rarefied OTU table was generated at 342 sequences per sample. Chao-1 and Shannon's H indices in fungal communities varied from 13 to 28 and from 1.904 to 2.598, respectively (Table 1). Both Chao-1 and Shannon's H showed the highest diversity in October, when chytrid cells showed the largest density (Fig. 2). Gephi network analysis is shown in Fig. 3. Results showed that the number of OTUs was very different in all months, and unique OTUs appeared in each month. The highest value of 11 unique OTUs was found in October, while the lowest value of 3 unique OTUs was found in November.

Temporal variations in zoosporic fungal communities
Temporal variations in zoosporic fungal communities were found during the study period. At the level of order, although the composition of the OTUs did not differ significantly throughout the study period, their sequence compositions showed temporal changes (Fig. 4). The sequence composition of Chytri diales, Rhizophydiales, Spizellomycetales and Monoblepharidales varied from 7.3 to 72.8%, from 13.0 to 56.2%, from 3.3 to 26.7% and from 0 to 19.8%, respectively, during the study period. Chytridiales dominated in May, June and September (> 65%) but not in other months (< 35%). Spizello mycetales were high contributors (> 21%) from October to December but not in other months (< 9%). Cladochytriales accounted for > 27% in November but were not detected from June to September. Only 3 known genera (Chytriomyces, Boothiomyces and Rhizo phydium) were detected throughout the study period. Six genera (Allomyces, Dendrochytridium, Rhi zoclosmatium, Gromochytrium, Gonapodya and Monoblepharis) were detected on only 1 occasion.
Three OTUs were matched with 3 chytrid species parasitic on phytoplankton: S. oculus as an obligate parasite of the green alga Staurastrum sp. (Van den Wyngaert et al. 2017), P. sphaericum as an obligate parasite of the diatom A. granulata (Seto & Degawa 2018) and Zygophlyctis asterionellae as an obligate parasite of the diatom A. formosa (Seto et al. 2020). The sequences of these 3 chytrids were always detected when sporangia of chytrids on each host phytoplankton species were found during microscopic examination (Fig. 6). However, their sequences were sometimes not detected even when chytrid sporangia

DISCUSSION
We detected diverse zoosporic fungi with the LSU region as a marker gene: a total of 101 OTUs in 3 zoosporic fungal phyla (Blastocladiomycota, Chytridiomycota and Cryptomycota) and 8 known species (Table 2). In our previous study based on the ITS region, many of the sequences were not assigned or were identified as 'unclassified fungi' (Song et al. 2018). Only 1 zoosporic fungal phylum (Chytridiomycota) with 35 OTUs was detected, and 1 parasitic chytrid (Staurastromyces oculus) was identified. Compared with the LSU region, there is less information on Chytridiomycota in the reference database of the ITS region for aquatic fungi (Wurzbacher et al. 2016).
Due to lack of chytrid data in public sequence databases, few OTUs were identified as known chytrids in the present study. The remaining OTUs were mostly assigned as unknown taxa, and therefore it is difficult to determine their ecology (Frenken et al. 2017). Eight OTUs, however, were identified as previously known species, and we were able to determine their ecological characteristics (i.e. their nutritional modes, whether they were parasites or saprotrophs, their hosts/substrates and habitats) ( Table 3).
Three OTUs were matched with 3 parasitic chytrid species of phytoplankton: S. oculus (Rhizophydiales, Staurastromycetaceae) infecting Staurastrum sp. (Van den Wyngaert et al. 2017), Pendulichytrium sphaericum (Chytridiales, Chytriomycetaceae) infecting Aula coseira granulata (Seto & Degawa 2018) and Zygophlyctis asterionellae infecting Asterionella formosa (Seto et al. 2020). The morphology of sporangia attached to Staurastrum dorsidentiferum was consistent with that of Staurastromyces oculus: it had an encysted zoospore in the isthmus region of the host, while sporangia attached on Staurastrum rotula were different from any stages of S. oculus (Van den Wyngaert et al. 2017). Different patterns in the abundance of sporangia on S. rotula and the relative abundance of the sequence reads of S. oculus (Fig. 6) also indicated that the chytrids infecting S. rotula must have been different species. The chytrid infecting S. dorsidentiferum may have been S. oculus, although there may have been other species infecting S. dorsidentiferum when sequence reads of S. oculus were not detected.
According to microscopic analysis, the morphology of sporangia attached to A. granulata was consistent with the morphology of P. sphaericum: the mature zoosporangium was spherical and stalked (Seto & Degawa 2018) (Fig. S3), and the patterns of read percentage and abundance of sporangia also matched (Fig. 6); thus, it is likely that the chytrid infecting A. granulata was P. sphaericum. The morphology of sporangia attached to A. formosa was consistent with that of Z. asterionellae, as the encysted zoospore was spherical and the developing thallus was obpyriform   (Fig. S3). Sequences of Z. asterionellae were always detected when chytrid sporangia on A. formosa were observed. The peak of Z.
asterionellae read percentage was a month earlier than the abundance of sporangia (Fig. 6). Since our DNA analysis targeted zoospores (2−20 μm), instead of sporangia, it is likely that the detected timing recorded a slight gap between zoospores and sporangia. Rhizophydium planktonicum is also a parasite of A. formosa and possesses similar morphology to Z. asterionellae ), but it was not detected with metabarcording analysis in this study. Some other parasitic chytrids might be included in the unknown OTUs. We found high infection rates in Closterium aciculare in July (Fig. 2) but could not identify any OTUs corresponding to parasitic chytrids on C. aciculare. Although chytrids on C. aciculare have been reported previously (Holfeld 1998, Fernández et al. 2011, no DNA sequence data are available for these chytrids. Among all the detected Chytridiomycota sequences in our metabarcoding analysis, 1 unknown chytrid OTU ('denovo 49959') might have been a parasite of C. aciculare, which showed relatively high abundance in July when the number of chytrid sporangia on C. aciculare was at its highest (Fig. 4). Phylogenetic analysis indicated that denovo 49959 may be a novel clade in the Rhizophydiales (Fig. 5). Further studies are necessary to clarify this, and to identify the parasitic chytrid of C. aciculare in Lake Biwa.
In the present study, the chytrid infection rates on 5 large dominant phytoplankton varied along with the abundances of their hosts. The temporal variation in hosts can be strongly influenced by the parasitic fungi, and especially by that of host-specific fungi (Holfeld 2000. In fact, chytrid infection has been shown to affect phytoplankton dynamics on S. dorsidentiferum in Lake Biwa . Chytrid infection-mediated effects on host phytoplankton dynamics may be plausible in this lake. Compared with the previous study , however, the prevalence of infection was low in our study. Since we only took monthly samples, we need more frequent sampling to identify the host− parasite relationships in the lake. The OTUs identified as Pateramyces corrientinensis, Globomyces pollinis-pini, Urceomyces sphaerocarpus, Rhizoclosmatium sp. and Pseudorhizidium endosporangiatum in the present study may be saprotrophic species, according to the literature (Table 3). They contributed a very small proportion of total chytrid reads (<10%), except in December (21%). P. corrientinensis, G. pollinis-pini and U. sphaerocarpus have been detected from pollen as substrates in aquatic habitats (Letcher et al. 2008). P. endosporangiatum has been isolated from pollen in Fig. 6. Comparison of known chytrids between metabarcoding analysis and microscopic counting results from 6 May to 12 December 2016 at Stn 3 in the north basin of Lake Biwa, Japan garden soil (Powell et al. 2013), and therefore in this study, it may have been transferred from terrestrial to aquatic habitats by wind action (Wurzbacher et al. 2010). Saprotrophic zoosporic fungi may depend strongly on substrate quality in the pelagic zone (i.e. carbon sources, sestonic C:N and N:P supply ratios from the watersheds) (Findlay et al. 2002, Güsewell & Gessner 2009), although we were unable to identify the factors influencing the dynamics of saprotrophic fungi in this study. Some OTUs showed the closest match with Neocallimastigomycota in this study (Table S3). This phylum is regarded as consisting of obligately anaerobic rumen fungi, which are found widely in the guts of ruminants (Gruninger et al. 2014). Both of our sequences matched 2 environmental sequences in GenBank, i.e. LLMB10_1 and LLME10_2 (accession numbers JN049554 and JN049555, respectively), from Lakes Tuscaloosa and Lurleen (Alabama, USA), respectively (Lefèvre et al. 2012), which were affiliated to Neocallimastigomycota. However, a BLAST search and preliminary phylogenetic analysis indicated that these sequences are not related to Neocallimastigomycota. LLMB10_1 was ca. 96% identical to the green algae Chlamydomonas spp. in the BLAST search. In our phylogenetic analysis (Fig. S4), LLME10_2 was not related to Neocallimastigomycota and showed a long branch in Chytridiomycota, indicating that the phylogenetic position of this sequence is uncertain. Wurzbacher et al. (2016) also found an OTU assigned to the Neocallimastigomycota sequence in Lake Stechlin (Germany), but concluded through phylogenetic analysis that the OTU was not related to Neocallimastigomycota but was distributed among other fungal lineages.
In our study, the sequences of the known chytrids were always detected by meta-barcoding analysis when they were observed by microscopic analysis. However, sometimes no sequences were detected even if the chytrid sporangia on phytoplankton were observed. One reason for this may be the nano-sized particle fraction (2−20 μm) that we targeted in metabarcoding analysis, but not in sporangia at tached to large phytoplankton species (> 20 μm). Metabarcod-  Table 2. Comparison of the detected zoosporic fungal phyla and orders using internal transcribed spacer (ITS) from our previous study (Song et al. 2018) and large subunit (LSU) regions by metabarcoding analysis. Samples were collected at Stn 3 in Lake Biwa, Japan. Plus (+) represents occurrence. Blank cells therefore indicate that these orders were not found. OTU: operational taxonomic unit  Table 3. List of identified chytrid species detected from 6 May to 12 December 2016 at Stn 3 in the north basin of Lake Biwa, Japan, including their ecology and hosts. -/-: not reported ing analysis was therefore unable to detect chytrid sequences when they mainly presented as sporangia on large phytoplankton. Although calcofluor white stains not only chytrids but also other fungi and protists with cellulose, we observed the morphology of sporangia carefully, and concluded that the detected sporangia were chytrids. The size threshold for chytrid infection has been uncertain, although chy trids tend to infect large hosts (Kagami et al. 2007a). In this study, we never found any chytrids that infected small phytoplankton, although it could be possible that many unexplained OTUs could imply miscounting those attached to small phytoplankton (≤ 20 μm). Another reason for this may be the primer bias, since we detected only a few fungi, and most of the sequences belonged to other eukaryotes. Further studies combining microscopic observations and metabarcoding analysis may reveal the existence of more diverse fungi in the water column of Lake Biwa.