Co-occurrences enhance our understanding of aquatic fungal metacommunity assembly and reveal potential host–parasite interactions

ABSTRACT Our knowledge of aquatic fungal communities, their assembly, distributions and ecological roles in marine ecosystems is scarce. Hence, we aimed to investigate fungal metacommunities of coastal habitats in a subarctic zone (northern Baltic Sea, Sweden). Using a novel joint species distribution model and network approach, we quantified the importance of biotic associations contributing to the assembly of mycoplankton, further, detected potential biotic interactions between fungi–algae pairs, respectively. Our long-read metabarcoding approach identified 493 fungal taxa, of which a dominant fraction (44.4%) was assigned as early-diverging fungi (i.e. Cryptomycota and Chytridiomycota). Alpha diversity of mycoplankton declined and community compositions changed along inlet–bay–offshore transects. The distributions of most fungi were rather influenced by environmental factors than by spatial drivers, and the influence of biotic associations was pronounced when environmental filtering was weak. We found great number of co-occurrences (120) among the dominant fungal groups, and the 25 associations between fungal and algal OTUs suggested potential host–parasite and/or saprotroph links, supporting a Cryptomycota-based mycoloop pathway. We emphasize that the contribution of biotic associations to mycoplankton assembly are important to consider in future studies as it helps to improve predictions of species distributions in aquatic ecosystems.


Introduction
Aquatic fungi comprise a diverse group of heterotrophic microorganisms, spanning a wide range of life strategies from saprotrophy through mutualism to parasitism (Nilsson et al. 2019). Their contribution and impact on the biogeochemical processes of the biosphere is crucial (Grossart and Rojas-Jimenez 2016, Gladfelter et al. 2019, Grossart et al. 2019. Molecular analysis of environmental samples using next-generation sequencing have unravelled an extremely high diversity of undescribed fungi (see e.g. Richards et al. 2015) that is often referred as the 'dark matter fungi', or DMF (Grossart et al. 2015). DMF are commonly found within the early-diverging lineages of the fungal tree of life, including members of the basal phyla of Cryptomycota and Chytridiomycota. These fungi have shown to be saprotrophs and obligate or facultative parasites, however, very limited knowledge is available about their actual ecological roles (Grossart et al. 2019). Saprotrophs participate in decomposition processes of dead organic matters, while parasitic fungi colonize living hosts to obtain nutrients for their life cycle (Donk and Bruning 1992). Parasitic fungi have a special significance because they can infect inedible phytoplankton species and facilitate energy transfer (Gleason et al. 2015) to zooplankton via their zoospores (Kagami et al. 2007, Agha et al. 2016, Garvetto et al. 2019), a mechanism known as mycoloop (Kagami et al. 2014). In turn, parasitic chytrid outbreaks can be mitigated and suppressed by grazers who feed on their zoospores, channelling nutrients across trophic levels (Kagami et al. 2011, Frenken et al. 2020. However, the field of aquatic chytrid biology is highly skewed toward culture-based studies, especially ignoring DMF, which prevents the achievement of a more complete, detailed understanding of this fungal group (Laundon and Cunliffe 2021). In marine ecosystems, our understanding of fungal host-parasite systems is even more limited, potentially due to the challenges in isolating and culturing the fungal partner (Gladfelter et al. 2019), and in determining their ecological roles which can vary under different circumstances (Grossart et al. 2015).
Aquatic fungi, and marine fungi in particular, have received much less attention in planktonic research compared to other planktonic microbes such as bacteria or protozoa . Despite that our knowledge on their diversity in marine ecosystems is greatly limited (Richards et al. 2015), a few studies have shown that distinct marine habitats harbour diverse and discrete fungal (mycoplankton) communities (Jeffries et al. 2016). Specifically, coastal habitats constitute a transitional zone between riverine and open ocean (offshore) sites. These coastal habitats function as sinks for terrestrial-sourced fungi (Hassett et al. 2019), and host higher proportions of the early-diverging fungal groups such as chytrids or Cryptomycota compared to oceans that are dominated by fungi belonging to Dikarya (Picard 2017, Hassett et al. 2019. Thus, coastal areas with high terrestrial influence should be hotspots for aquatic fungi. The high fungal diversity and the elevated proportion of chytrids in the Baltic Sea may support this hypothesis (Hassett et al. 2019, Rojas-Jimenez et al. 2019. Nevertheless, we have limited knowledge on how compositional shifts from freshwater inlets towards offshore sites differentiate mycoplankton communities, and what processes (e.g. community assembly mechanisms) affect their distributions.
Metacommunity concept provides a framework to reveal important mechanisms that simultaneously shape (and maintain) biodiversity variation at regional scale (Leibold et al. 2004). These mechanisms include, for instance, selection by the environment (environmental filtering), dispersal-related processes (e.g. dispersal limitation, habitat connectivity) and stochasticity (ecological drift). To our knowledge, only a single study (Yang et al. 2021) aimed to estimate these processes in mycoplankton communities and found that their influences significantly differ along a riversea transect, having more stochasticity in the coastal and offshore sites, while environmental selection dominates in upstream sections. These community assembly processes can have distinct influences on each species and in each habitat (Leibold et al. 2022). Furthermore, the influence of biotic interactions on metacommunity structure, which was lacking in previous concepts, can be important and, thus, should be acknowledged (Leibold et al. 2022). Joint species distribution models (JSDMs) have the advantage to provide such important insights into the drivers of variation in species distributions, and with this, estimate how the contributions of space, environment, and biotic interactions, driving metacommunity assembly, differ among sites and species (Ovaskainen et al. 2019, Poggiato et al. 2021, Leibold et al. 2022). In addition to distribution models, network approaches have been applied successfully to reveal detailed biotic associations. Such methods were proven to predict host-parasite interactions in virus-and chytridfocused studies (Rojas-Jimenez et al. 2017, Kilias et al. 2020, Meng et al. 2021, Ilicic et al. 2022) as well, even if network inference needs to be taken cautiously. Nevertheless, application of network approaches can help us to extend current knowledge on the relationships of fungi with other microorganisms (Rojas-Jimenez et al. 2017).
In this present study, we aimed to investigate fungal (mycoplankton) metacommunities of coastal habitats in a subarctic zone (Gulf of Bothnia, northern Baltic Sea, Sweden). We applied long-read metabarcoding approach using a Nanopore Min-ION sequencing platform which has been shown to provide accurate and efficient data in the characterization of aquatic environmental DNA (Davidov et al. 2020). Using recently developed novel joint species distribution model and network approach, we aimed to quantify the importance of biotic associations contributing to fungal metacommunity assembly, further, detect potential biotic (i.e. parasitic) interactions between fungi-algae pairs, respectively. We hypothesized compositional shifts of fungal metacommunities along the freshwater inlet-bay-offshore transects, with high fraction of early-diverging fungi, especially chytrids, in the bays. Since the results of the joint species distribution model would allow us to infer community assembly processes (Leibold et al. 2022), we assumed that mycoplankton communities of bays should be influenced by a high degree of ecological drift rather than environmental selection. Furthermore, we expected that cooccurrence network reveal putative chytrid-algal interactions, as well as new links with Cryptomycota fungi.

Sample collection
Samplings were conducted by following the course of four freshwater inlets towards the brackish water of the Gulf of Bothnia, northern Baltic Sea. Four surface water samples were collected (0.5 m depth) once in a month in these four bays during the summer period (June to September) in 2018 (Eriksson et al. 2022). In addition, samples from two offshore sites were also collected in each month (Fig. 1). All water samples were collected in sterile sampling bottles, then transported to the laboratory and stored in the dark at 4 • C. Physicochemical properties of samples were measured (Eriksson et al. 2022) and included in this study. Briefly, temperature, pH, and salinity were measured in situ with a WTW ProfiLine Cond 3110 portable device, and subsamples were used to measure environmental variables including dissolved organic carbon (DOC), total dissolved nutrients (TDP, TDN), dissolved inorganic nutrients (e.g. PO 4 3− , NH 4 + , NO 2 − , NO 3 − , SiO 2 ), and humic substances. Environmental data for the inlet (freshwater) samples of each bay is lacking. For the microbial community analysis, we aimed to filter 500 ml subsamples through sterilized 0.2 μm pore size membrane filters (Supor, 47 mm, Pall Corporation) the same day as the water collection and kept the filters at -80 • C until further processing.

DNA extraction and sequencing
The DNA was extracted from the filters using the DNeasy Pow-erWater Kit (Qiagen) according to the kit protocol with some modifications. Namely, the samples were treated with an additional heating step (applying horizontal water bath for 30 min at 65 • C) and with a bead-beating step in each direction of 20 Hz, 3 × 3 min with a TissueLyser II (Qiagen) in order to aid the lysis of microorganisms including fungal cell walls. DNA extracts were quantified with a Qubit 4 fluorometer (Thermo Fisher Scientific).
To perform metabarcoding, we used the NS1short/RCA95m primer pair in order to amplify the majority of the ribosomal operon of the ribosomal tandem repeat (18S-ITS1-5.8S-ITS2-28S rDNA operon) ). The PCR was performed, similarly as in Wurzbacher et al. (2019), but see details in Supplementary material. Amplification and its fragments length (∼4-6 kbp) were confirmed by gel electrophoresis resulting in 61 final samples (out of 72). The barcoded PCR products were purified with 0.8× of AMPure magnetic beads (Beckmann) following the manufacturer's protocol. Thereafter, the purified PCR products were quantified using the Qubit 1× HS Assay Kit (ThermoFisher Scientific), pooled in equimolar amounts. This final pool was concentrated with 2.5 × of AMPure magnetic beads in 50 μL nuclease-free water (ThermoFisher), and again quantified (Qubit 1× HS Assay Kit).
1 μg of library was used for the ONT library preparation using the 1D sequencing (SQK-LSK109; Oxford Nanopore Technologies). Sequencing was done on a MinION Mk1C instrument (ONT) operated with a Spot-ON Flow Cell (R9.4.1 chemistry). Real-time basecalling was executed using the High-accuracy basecalling (HAC) mode using the MinKNOW software (v21.05.12). In the end, 1.42 M reads were yielded with N50 read length of 4.4 kb and Q > 9.

Sequence data processing
Quality reads were demultiplexed and barcoded primers were trimmed with MiniBar (Krehenwinkel et al. 2019), and then filtered by length (3-7 kb) with NanoFilt (v2.8.0) (De Coster et al. 2018). These filtered, quality reads were then processed in the recently developed NGSpeciesID pipeline (v0.1.2.2) which wraps a set of tools to generate clusters and form polished consensus sequences for each cluster (Sahlin et al. 2021). It used the isONclust algorithm (Sahlin and Medvedev 2020) for read clustering (-mapped_threshold 0.8 -aligned_threshold 0.5) which accounts for variable error rates within reads. Draft consensus sequences were formed for each cluster containing at least five reads with spoa (v4.0.7) using maximum 500 sequences (-max_seqs_for_consensus 500), then reverse-complement clusters were merged using Parasail (v1.2.4). Finally, polishing of the consensus sequences (input reads are mapped back to the consensus sequence and basepair errors are corrected) were done with Racon (v1.4.20;Vaser et al. 2017) using two iteration steps. Since polished consensus sequences are the final output of the NGSpeciesID pipeline, we formatted the output files using custom scripts in R (R Development Core Team 2016) to have an appropriate input file for Mothur (v1.46.1; Schloss et al. 2009) to assign the corresponding sequence count to each OTU. Sample coverage was assessed with the 'iNext' R package (Hsieh et al. 2016), and found that, with the exception of four samples, community composition was sufficiently covered (Fig.  S1).
Taxonomic classification of the 512 final consensus sequences was primarily done by local BLAST search against the NCBI nucleotide database (nt) (downloaded on 23 October 2021) using BLAST+ (v2.11.0+), limiting the BLASTn search for Fungi (taxids: 4751) and keeping hits with at least 80% identity. Results of the BLASTn search was processed with phyloR (https://github.com/c parsania/phyloR) to keep top hits (Supplementary Table S2) and to assign taxonomy levels. To support taxonomy, rRNA genes (SSU and LSU) and the entire internal transcribed spacer (ITS) were extracted using ITSx (Bengtsson-Palme et al. 2013), and used in BLASTn search to assign taxonomy against secondary reference databases: SILVA Reference SSU and LSU databases (Quast et al. 2013) (release 138.1), and the UNITE+INSD database (Abarenkov et al. 2021) (version 10.05.2021) in case of ITS (see , Table S3). Due to the prevailing classification challenges and conflicts among these secondary databases (Heeger et al. 2018), we only used this secondary classification step to validate fungal sequences as follows. The OTUs that had not matched with any fungal sequences in the secondary reference databases were discarded. In cases when OTU was classified as fungi by only one of the secondary databases, we manually choose the classification with the highest bit score. OTUs without matches below Kingdom level based on the BLASTn search against NCBI nt database were classified as 'likely fungi' if all secondary classifications would match fungal sequences. In the end, 493 OTUs were identified as fungi (with an average alignment length of 2130 bp and 93.2% identity), thus, were kept for downstream analyses (Table S4). Potential chimeric consensus sequences identified by ITSx and their corresponding OTUs from the OTU table were filtered out.
Raw reads were deposited to NCBI SRA database under the accession number PRJNA849821.

Data analyses
Venn diagram was used to visualize the number of shared fungal taxa across sampling sites, furthermore, diversity analyses (alphadiversity and beta-diversity based on Bray-Curtis distance) were performed using the 'microeco' R package (v.0.6.5) (Liu et al. 2021) and the results (i.e. nonmetric multidimensional scaling-NMDS) were plotted using 'ggplot' package (Wickham 2009). Difference in alpha diversity across bays, their inlets and offshore sites were tested with ANOVA followed by Duncan's test (P < 0.05) as a posthoc test. To test compositional differences between samples, permutational multivariate analysis of variance (PERMANOVA) with 999 permutations was performed using the function adonis2 in 'vegan' R package (Oksanen et al. 2016). Distance-based redundancy analysis (dbRDA) was also performed on the bay samples to assess the influence of environmental variables or sampling time (i.e. Day). To statistically test their influences, Mantel test (Spearman correlation with 999 permutation) were computed on community dissimilarity based on Bray-Curtis distance.

Inference of the internal structure of metacommunities using joint species distribution modelling
To predict metacommunity structure as a whole, a recently developed scalable joint species distribution model (sjSDM) approach (Pichler and Hartig 2021) was applied which use Monte Carlo integration of the JSDM likelihood together with elastic net regularization on all model components. Due to the lack of measured environmental data for freshwater inlets and offshore sites, we used only samples from the four bays (n = 40, no. of OTUs = 442). Spatial eigenvectors were generated from the GPS coordinates to account for spatial autocorrelation (using the generateS-patialEV function) and measured environmental variables were z-transformed prior the analysis. The regularization for all covariances and coefficients (following general suggestions from the developers) was tuned over 40 random steps with leave-oneout cross validation (LOOCV), 150 iterations and learning rate of 0.01 using the sjSDM_cv function of the sjSDM R package (v1.0.1) (Pichler and Hartig 2021). Thereafter, the best regularization parameters were used to fit a multivariate probit model using the sjSDM.tune function. For this analysis, the Python library 'PyTorch' (Paszke et al. 2019) was run from within R thanks to the 'reticulate' R package (Allaire et al. 2018). The model is available as R Data file in Open Science Framework (OSF) (http://osf.io/764mu).
Following the framework by Leibold et al. (2022), our model was then used to estimate how the contributions of environment, space and biotic associations shaping metacommunity assembly vary among sites and taxa. This quantitative approach allowed us to identify how different sites and taxa contribute to the overall metacommunity structure, taking into account the fact that a set of sites or taxa are not necessarily equally influenced by the interplay of environmental filtering (abiotic selection), dispersal, biotic interactions and ecological drift (Leibold et al. 2022).

Network analysis
With the aim to further investigate biotic interactions and, in particular, potential parasitic-host associations, co-occurrence network was constructed. The 18S V6-V8 rRNA gene sequences from the study of Eriksson et al. (2022) and the subset of their dataset (deposited in OSF; http://osf.io/764mu) were used to assess potential host-parasitic and/or saprotrophic associations in our bay samples (n = 40). Absolute read counts were used for the network analysis because previous studies have shown that relative abundance data suffer from apparent correlations which lower specificity of association networks (Berry andWidder 2014, Meng et al. 2021). Further, since low number of sites are susceptible to false positive correlations (Berry and Widder 2014), we decided to run one global network analysis without subsetting our dataset into the three time periods or to individual bays. Nevertheless, in this analysis, we used only the most dominant fungal groups (Cryptomycota and chytrids) that are prone to have parasitic lifestyle, but also included the 'likely fungi' group to assess their putative cross-kingdom associations. To account for variations in sequencing depth between the long-read fungal and 18S-based algal datasets, we decided to use FlashWeave (v0.18.1) (Tackmann et al. 2019) which applies clr transformation to handle compositionality with its adaptive pseudo-counts (clr-adapt) approach. Julia (v1.6.4) was used for constructing association network with FlashWeave package from within the 'microeco' R package (Liu et al. 2021).
Prediction of ecological interactions (P < 0.01) between fungi and algae was performed using FlashWeave-S ('sensitive' mode) with default settings. For enhanced reliability, associations were computed only when an OTU was present more than 20 times (automatically determined by the software). This approach also deals with indirect associations (as a result of shared environmental niches) by the incorporation of environmental metadata. The integration of 15 meta-variables (MVs such as 'Day' as sampling time, 'Bay' as identities, and measured physicochemical variables) was done in order to remove potential indirect associations (i.e. as a result of shared niche preference). The constructed cross-kingdom network was visualized with Gephi (v0.9.2) using associations' correlation weight > |0.4| and the ForceAtlas2 layout (Jacomy et al. 2014). For clarity MV nodes were moved to the side.
Fungal community compositions exhibited strong compositional shift over time (from June to September, F = 11.36, P = 0.001), although the PERMANOVA test also confirmed spatial variability (F = 1.79, P = 0.025) which was mainly pronounced by the separation of bay samples from their respective freshwater inlets (Fig. 4a). The significant interaction between sampling sites and time (F = 1.96, P = 0.005) further showed that temporal variation differed between sites. When the fungal metacommunities of the four bays were analysed alone (using dbRDA) to investigate the influence of the measured environmental variables (Fig. 4b, Table S6), temperature, pH, humic substances, salinity, DOC, NH 4 + , TDN, and, to a lesser degree, NO 3 − and chlorophyll-a concentrations were important (all P < 0.05) in shaping the beta diversity of fungal metacommunities, besides sampling time (i.e. Day) (see the results of Mantel test in Table S6). The three distinct clusters are clearly apparent from the dbRDA (Fig. 4B) and accord with sampling time rather than with space (i.e. bays) when environmental parameters were considered. These temporal dynamics in fungal metacommunities can be attributed to several trends observed in the relative abundances of different taxonomic groups (see Fig. S3). Specifically, Cryptomycota (e.g. Paramicrosporidium sp. and mainly unidentified Cryptomycota (73.1%)), Ascomycota (members in the genera of Articulospora, Cladosporium, Penicillium), Basidiomycota (Vishniacozyma sp., Rhodotorula sp.) and Microsporidia (Mitosporidium sp.) showed elevated abundances towards late summer (August and September), while Chytridiomycota had its peak in July (dominated by OTUs assigned as Betamyces spp. and Chytridium polysiphoniae) and decreased similarly as 'likely fungi' by the end of the sampling campaign (September). In the remaining fungal groups, there was no clear, detectable trend.

Internal structure of aquatic fungal metacommunities
We used scalable joint species distribution model to estimate the importance of environment, space and biotic interactions in driving metacommunity assembly among sites and fungal taxa of the bays. Our model revealed that the distribution of fungal taxa was mainly driven by environmental factors and to a lesser extent by spatial and biotic associations (Fig. 5, Table S7). Sites within each bay differed greatly in how their community compositions were attributed to the three components, and were explained by R 2 values ranging from 0.4% to 1.6% (Fig. 5, left panel; Fig. S4). Importance of biotic interactions (or co-distributions) was enhanced in several samples regardless of the sampling time or bay (Fig. S4, Table S7).
Taxa distributions were weakly predicted by the model (e.g. low R 2 values) and were stronger influenced by biotic associations as the effect of the environment decreased (Fig. 5, right panel). The model, however, could predict taxa distribution to a greater extent when environment had strong influence. In contrast, lower predictive power occurred when taxa were almost exclusively influenced by space or biotic associations. OTUs within each dom-inant phylum were distinctly affected by the three components (Fig. S5), hence, the distribution of each taxon was determined by a unique combination of environment, space and biotic factors. Phyla that are represented by a very few members (i.e. Blastocladiomycota, Microsporidia, Mucoromycota and Zoopagomycota), and hence considered as rare, were mainly influenced by environmental factors, while a Sanchytriomycota occurrence was mainly affected by spatial factors. It is, however, important to mention that temporal factor (i.e. Day), despite its strong influence as presented in the above-mentioned multivariate analyses, was not included in the model in order to investigate the pure effects of measured environmental variables on taxa distributions.

Co-occurrence network
We generated one cross-kingdom co-occurrence network on all bay samples to predict ecological interactions between specific fungal groups (Cryptomycota, chytrids and 'likely fungi') and algae (Fig. 6, S6). In total, this final network was composed of 288 nodes (276 OTUs) and 364 associations (edges). Detailed list of all interactions can be found in Table S8. Among the three fungal groups (Cryptomycota, Chytridiomycota and 'likely fungi') 120 interactions were detected, of which 92.5% were positive. The strongest associations (e.g. corr. weight > +0.9) occurred between OTUs within the same taxonomic group. Between Cryptomycota and Chytridiomycota six associations (one negative and five positive) were found, and these chytrid OTUs had no positive link to any algal OTUs. To reveal potential parasitism (corr. weight > +0.4) between kingdoms, we found 25 interactions. Most (13) were identified between chytrids and algae (five Ochrophyta, four Chlorophyta, three Dinoflagellata, and only one Cryptophyta).    Measured environmental parameters and/or sampling time can lead to spurious correlations between OTUs that share similar niche optima and thus associated with the same meta-variable (MV). The inclusion of MVs in our network analyses showed that several fungal (8) and algal (18) OTUs were mostly associated with water temperature and chemical properties (e.g. nitrite and humic substances), respectively. Sampling time (Day), interestingly, only associated with three algal taxa: Nannochloris sp. (Chlorophyta) showed positive association (e.g. positive trend in their abundances over the sampling course), while Diatoma tenuis and Uroglena americana had both negative associations. There was a Cryptomycota OTU showing positive association with sampling sites (Bay) and this OTU was also highly influenced (40.1%) by spatial effects as revealed by sjSDM.

Composition and diversity of coastal mycoplankton
Our long-read metabarcoding study identified 493 fungal taxa along four transects, encompassing freshwater inlets up to their corresponding bays and beyond to offshore habitats in the Gulf of Bothnia, Sweden. Generally, alpha diversity of mycoplankton declined and community compositions changed along these four transects (from freshwater to offshore sites). This resonates with a declining trend found by Yang et al. (2021) in a much greater spatial scale along the Elbe River down to its estuary (North Sea). The general compositional shift over sampling sites, however, differed from the one found in the study of Yang et al. (2021). Namely, the authors found greater dominance of chytrids in samples with greater freshwater influence and the dominance of Ascomycota and Basidiomycota in open sea (offshore) environments, while our results suggest the opposite. This disparity may arise due to a greater salinity range their study covered, increased anthropogenic activities around that study area (i.e. influence of the city of Hamburg), or due to the lack of the inclusion of temporal aspect in their study. All these might affect mycoplankton communities in space, creating distinct community structures and disparities among studies.
A great fraction of the mycoplankton (44.8%) was assigned as early diverging fungi (i.e. Cryptomycota and Chytridiomycota) which aligns with previous studies that reported the dominance of these fungal phyla in both freshwater and marine environments (Comeau et al. 2016, Hassett andGradinger 2016). This indicates that early-diverging fungal phyla dominate not only sheltered aquatic ecosystems, as previously suggested (see e.g. Rojas-Jimenez et al. 2017), but also in habitats which are greatly affected by terrestrial input and continuous mixing events. As we hypothesized, chytrids showed elevated richness in bays and offshore sites, and dominated in mid-summer (i.e. July). Dominance of chytrid-like sequences was also found in a field survey, covering a great selection of marine (Arctic) and freshwater (temperate biomes) environments (Comeau et al. 2016). Their results, however, showed much less frequency and relative abundance of Cryptomycota, in contrast to our findings. Picard (2017) has found greater dominance of Ascomycota in mycoplankton communities in coastal North Carolina, and chytrids were more dominant fraction of the fungal communities in sediments, suggesting the importance of dormancy (Velasco-González et al. 2020). Previous studies targeting aquatic microfungi are scarce, but the findings so far indicate that the proportion of different fungal groups can be greatly dependent on the spatial scale and the environment of interest.
In spite of the high diversity of 'dark matter fungi' in aquatic environments (Grossart et al. 2015), we clearly lack detailed taxonomic information for most Cryptomycota OTUs (88.2%), even when single marker genes are targeted. Although, their ecological importance is most likely relevant. This may not be surprising, giving the fact that newer clades and taxa are continuously discovered and described, thus, the taxonomy of such earlydiverging fungal lineage is not completely resolved. Among Cryptomycota, an intranuclear parasite of amoebae (Paramicrosporidium saccamoebae) occurred in most freshwater inlets, further, with an elevated abundance during August and September in the bays. Even though freshwater inlets could nearly continuously support dispersal of parasitic Cryptomycota into coastal and offshore habitats, its establishment is most likely determined by other factors (e.g. environmental conditions and/or hosts presences). Interestingly, two occurring Rozella spp. (commonly described as hyperparasites) did not end up in the co-occurrence network, suggesting the lack of hyper-parasitism and/or algae parasitism of these taxa in the observed bays.
Betamyces spp. and Chytridium polysiphoniae were dominant members of chytrids in July and in the bays, particularly, and have been reported as parasites on several phytoplankton species (Christaki et al. 2017) and seaweed in estuarine and marine habitats (Gleason et al. 2011). This, in line with our hypothesis, further supports that coastal habitats provide an optimal environment for chytrids (Hassett et al. 2019) whereas the mixing of riverine and saline waters promotes the growth of phytoplankton and source aquatic fungi with parasitic lifestyle. The peak relative abundances of chytrids in July may indicate the presence of parasitic taxa which usually emerge at elevated algae biomass, channelling nutrients to higher trophic levels via the mycoloop (Grami et al. 2011, Gerphagnon et al. 2015, Frenken et al. 2017. Numerous members of Ascomycota and Basidiomycota (e.g. Cladosporium and Blumeria) found in this study are known as fungi with terrestrial origin. These taxa potentially represent nutrient sources for saprotrophic aquatic fungi (Magyar et al. 2016). However, more and more evidences suggest that a fraction of these fungi display a truly amphibious ability, hence, their presence as mere metabolically inactive spores or relictual DNA is questioned , Koivusaari et al. 2019. Phylum Ascomycota was dominated by Articulospora atra, which is a commonly known aquatic hyphomycetes, identified in numerous freshwater ecosystems, previously. To our knowledge, this is the first record showing that it can establish populations in brackish ecosystem, likely dispersed from the freshwater inlets into the bays.

Metacommunity dynamics
Besides the major impact of salinity in shaping fungal communities (Rojas-Jimenez et al. 2019, Ilicic et al. 2022, our results highlight the relevance of other abiotic factors (i.e. temperature, pH, humic substances, etc.) as well as the influence of time in affecting community dynamic of mycoplankton. This emphasizes that snapshot studies easily can miss important aspects of community dynamics, which are clearly influenced by the simultaneous influence of space and time. In general, the importance of assembly processes shaping aquatic fungal communities is largely unknown. To our knowledge, only Yang et al. (2021) aimed to estimate community assembly processes and found that selection was more important in less saline environments (freshwater inlets) together with dispersal limitation, while communities of the coastal and offshore sites were mainly assembled by mass (dispersal) effect and ecological drift. By inferring the internal structure of metacommunity using modelling approach, our results suggest that environmental filtering had strong influence on the observed bay-inhabiting microfungi, and dispersal was not limited as evidenced by much lower importance of spatial effects than environment. In extreme cases when taxa were highly influenced by space and/or biotic associations, their distributions were predicted very weakly by the model, as the low R 2 values suggest. This supports Leibold et al.'s (2022) simulation results that suggested greater predictive power for species with narrower (more specialized) environmental niche than those that are dispersal limited. The generally low predictive power indicates that taxa-level distribution of mycoplankton is greatly influenced by ecological drift (stochasticity), hence, more challenging to predict compared to community-level distributions. Although, the low predictive power may also be the result of several other factors (i.e. trophic interactions, unmeasured environmental parameters) that were not included in our model.
Our overall findings suggest that both environmental and biotic factors can easily enhance processes (i.e. selection driven by abiotic environmental factors and/or competition) that contribute to the internal structure of metacommunities, determining the distribution of microfungi. The inclusion of biotic associations highlighted that their influence can strengthen when the abiotic and spatial drivers are not dominating. It is important to note that the biotic component of the model (i.e. biotic associations) might mask processes resulting from ecological drift (stochasticity), while spatial patterning might include unmeasured but spatially autocorrelated environmental variables (Blanchet et al. 2020, Wilkinson et al. 2021. In field studies, therefore, the presence of complex trophic interactions (i.e. viruses, grazers), unmeasured environmental variables may all represent additional challenges in teasing apart processes that structure microbial communities, and thus, cause biased detection of the effect of environment-space-biotic components in this particular joint species distribution model applied here. Still, by taking into account interactions (or more precisely covariances) between species, joint species distribution models can represent a useful tool for researchers who aim to investigate distributional patterns.

Putative biotic interactions
Our network analysis revealed a diverse picture, with a high number of co-occurrences among three fungal groups (Cryptomycota, Chytridiomycota, and 'likely fungi') which supports our sjSDM result that estimated in average 18% contribution of biotic associations among these groups. This highlights that the co-occurrences identified by network approaches can indicate direct pairwise interactions, although it might also suggest co-distributions via a complex interplay between true biotic associations, unmeasured environmental factors and, naturally, ecological drift (stochasticity) (see e.g. Wilkinson et al. 2021).
Twenty-five positive associations were found between fungi and algae. Interestingly, one of the most abundant chytrids (Betamyces spp.), which are recognized as parasites on numerous phytoplankton species (Christaki et al. 2017), showed strong associations with only an unassigned Dinophyceae taxon. Most interactions between chytrids and algae involved multiple partners which suggests that these potential parasitic interactions work in multifaceted ways to impact host distribution and abundance.
However, it is important to note that these potential pathogens might turn into saprotrophs upon host decay. In that scenario, our prediction of the ecological role of the associated partners can be greatly challenged (Egan and Gardiner 2016).
Previous studies suggested the importance of Cryptomycota (e.g. Rozella) species in the regulation of population size of parasitic zoosporic fungi in lakes (Gleason et al. 2012). We found, in contrast to our initial assumption, no sign of Cryptomycota OTU which may have been acted as hyper-parasites of parasitic chytrids (i.e. chytrid showing strong associations with algal OTUs). Although some strong links were found between Cryptomycota OTUs and Chytridiomycota OTUs, these chytrids were hypothetically saprotrophs as they lacked clear positive relationship with any algae. Hence, we cannot rule out the possibility that these Cryptomycota-Chytridiomycota links might indicate mere co-distributions or other biotic interactions. Nevertheless, detected associations between members of Cryptomycota (unassigned) and several algal phyla (i.e. Cryptophyta, Dinoflagellata) suggest that this fungal group might have a relevant role in the regulation of algal populations which are common in coastal habitats and in the Baltic Sea (Bazin et al. 2014, Enberg et al. 2018, without being targets for chytrids (Gleason et al. 2015). This, in turn, emphasizes that Cryptomycota may have as important role as chytrids in aquatic food webs, and reinforce the existence of other than chytrid-based mycoloops in aquatic ecosystems (Kagami et al. 2014, Priest et al. 2021).

Conclusions and future perspectives
This study enhances our knowledge of fungal diversity in coastal marine habitats, elucidate their spatiotemporal variation, and presents biotic (i.e. parasitic) interactions with planktonic algae. We believe that the joint application of distribution models and network approaches represents the advantage to infer a more detailed picture of metacommunity assembly with the inclusion of species covariation (attributed to biotic interactions), and can be used to support conclusions drawn from network results (e.g. presence of biotic interactions between Cryptomycota and chytrids).
Previous studies investigating fungal interactions have mainly been restricted to studies in laboratory settings. Whilst those studies deepen our knowledge in biotic interactions and their mechanisms, they barely provide information about how these interactions actually impact the distribution of mycoplankton communities in nature. Our findings emphasize that the contribution of biotic associations to fungal metacommunity assembly are important to consider in future studies as it helps us to improve predictions of species distributions in aquatic ecosystems. Identifying biotic relationships that affect the distributions of members of mycoplankton could be useful for plankton ecology, through habitat management to promote species which control algal blooms and facilitate nutrient transfer to upper trophic levels.
Studies targeting aquatic fungi, with a special focus on 'dark matter fungi', have a great potential. Therefore, future works should extend chytrid-centered studies on other fungal groups (e.g. Cryptomycota) and consider studying their ecological roles and community assembly in a wider selection of aquatic ecosystems. We also imagine future field studies as being the basis of laboratory experiments wherein putative biotic interactions, moreover, the ecological roles of mycoplankton members could be further assessed using more directed (e.g. single-cell genomics and meta-omics) approaches (Laundon and Cunliffe 2021). Taken together, we emphasize the need to go a step beyond culture-based studies and approach aquatic fungi from a (meta)community level perspective, too, in future studies.

Author contributions
AA, UCG, JW conceived the study and obtained funding for this research. KE performed sample processing and provided samples for this study. MV processed samples further (including sequencing), developed the research questions, analysed the data and drafted the manuscript. All authors contributed substantially to writing and revisions.