Non-random patterns of chytrid infections on phytoplankton host cells: mathematical and chemical ecology approaches

Host−parasite interactions between phytoplankton and fungi (chytrids) are key processes in aquatic ecosystems. However, individual-level heterogeneity in these interactions remains unexplored, although its importance in predicting the spread of diseases has been demonstrated in epidemiology. In this study, we experimentally tested whether individual-level heterogeneity could be a good indicator of phytoplankton−chytrid interactions, using a freshwater green alga Staurastrum sp., the diatoms Ulnaria sp. and Fragilaria crotonensis, and chytrid fungi. The number of attached fungi per host cell showed a non-random clumped parasite distribution on Ulnaria sp. and F. crotonensis, but a random Poisson distribution on Staurastrum sp. To explore the potential mechanisms of these patterns, we developed a mathematical model describing sequential encounters between chytrid zoospores and host cells. The statistical fits of the model explained the parasite distributions for Ulnaria sp. and F. crotonensis well, indicating that the clumped parasite distributions may result from an infection rate, increasing with the number of infections that already occurred on each host cell. Simultaneous analysis of volatile organic compounds (VOCs) from uninfected and infected host populations revealed that, among 13 VOCs detected, 6 components characterized the differences in VOC compositions between species and infection status. In particular, the level of beta-ionone, potentially acting against fungal activities, was significantly reduced in the presence of chytrid infection of Staurastrum sp. These VOCs are targets for future studies, which potentially act as chemical signals influencing chytrid zoospores’ behaviors. The combination of mathematical and chemical analyses represents a promising approach to better understand the individual-level processes of phytoplankton−chytrid interactions.


INTRODUCTION
Interactions between phytoplankton and fungal parasites such as chytrids are key processes in aquatic food webs and ecosystems, potentially driving phytoplankton population and community dynamics (Miki et al. 2011, Frenken et al. 2017, Song et al. 2021 ) and the flow of energy from phyto-plankton to zooplankton via zooplankton grazing on chytrid zoospores, the so-called 'mycoloop' (Kagami et al. 2007(Kagami et al. , 2014. Past studies have estimated the prevalence of infection in natural ecosystems (Ibelings et al. 2011, Kagami et al. 2012, experimentally quantified the interaction strength (Kagami et al. 2004, Frenken et al. 2020, and developed mathematical models (Bruning 1991a, Miki et al. 2011, Gerla et al. 2013, Almocera et al. 2018, Frenken et al. 2020. These studies demonstrate that the mycoloop represents one of the major trophic pathways in aquatic systems. However, most studies do not focus on individuallevel heterogeneity in the interactions between host cells and parasitic chytrids (except Bruning 1991b, Van den Wyngaert et al. 2014), but focus on population-averaged parameters only (e.g. the prevalence of infection, transmission, and of zooplankton grazing rates on chytrid zoospores). This is based on an implicit assumption that the encounter between host cells and chytrid zoospores occurs at random and that the dynamics of host−parasite interactions can be reasonably understood with the population-averaged parameters. Even with the existence of chemotactic behavior of chytrid zoospores (reviewed by Gleason & Lilje 2009) and species-specific or strainspecific attachment of zoospores (Kagami et al. 2007), its effect on population dynamics is understood as an increase in the encounter rate (see Lambert et al. 2019 for bacterial chemotaxis) and host-specific encounter rates. When chemotaxis is explicitly modeled as directed movement of zoospores depending on the local spatial distribution of host cells, it shows the importance of this behavior to explain parasite infection success, especially at low host density (Van den Wyngaert et al. 2014). However, individual host heterogeneity and impact on population dynamics remains less explored. Contrary to phytoplankton− chytrid interactions, human epidemiology and general parasitology consider the non-random interactions between the host and infective agents and nonrandom individual-level heterogeneity in hosts as major determinants of the spread and persistence of infectious diseases (Paull et al. 2012). Non-random individual-level heterogeneity has been studied in sexually transmitted diseases (STDs; Woolhouse et al. 1997), severe acute respiratory syndrome (SARS; Lloyd-Smith et al. 2005), and even computer viruses (Vazquez et al. 2007). In such diseases, the transmission rate (i.e. the potential of infected hosts to cause secondary infections) shows large variations between host individuals and there is a small fraction of host populations with a very high transmission potential (superspreaders). Then, even with a population-averaged low transmission rate, a limited number of superspreaders contribute to the spread and persistence of the infection in the host population. This is because such superspreaders disproportionally contribute to new infections and increase the basic reproduction number (R 0 ), which is the number of secondary infections expected by a single initial infection (Paull et al. 2012). This can also be seen in the current global COVID-19 pandemic (Lee et al. 2020).
In addition, non-randomness of host searching behaviors, which distinguish the state and quality of hosts, have been well studied in plant−herbivorous insect interactions and herbivore−natural enemy insect interactions in terrestrial food chains (Arimura et al. 2009). Theoretical approaches have long focused on non-random behaviors in resource−consumer interactions, which are demonstrated to be the basis of any functional response (Cosner et al. 1999), the determinant of population stability (Anderson & May 1978, Hollingsworth et al. 2015, and the key of their evolutionary consequences (May & Nowak 1995, Mideo 2009, Yoneya & Miki 2014. These past studies in other disciplines strongly indicate that non-random individual-level heterogeneity in phytoplankton− chytrid interactions -if existent -must be a potentially critical factor of both phytoplankton and chytrid population dynamics.
In this study, we therefore aim to demonstrate the existence of non-random interactions between phytoplankton hosts and parasitic chytrids. Then, we also aim to explore the potential mechanisms of such nonrandom interactions and to develop a quantitative method to describe their non-randomness. Specifically, we address the following 3 questions.
First, can we practically distinguish random and non-random interactions between host and chytrid by simple microscopic observations? Our strategy is to distinguish these interactions by monitoring the general infection patterns, instead of directly observing the infection processes. This is because it is not easy to directly observe all phases of the infection processes (e.g. the attraction of zoospores to a host, the maturation of sporangia, and release of newly produced zoospores) at the individual level. A key idea is that the existence of individual-level heterogeneity of infection patterns (e.g. the number of infections) does not immediately imply the non-randomness of infection processes. A simple solution is to count the number of attached fungi (attached zoospores, sporangia, and empty sporangia) per host cell, derive its frequency distribution (hereinafter, infection number distribu-tion) in a population, and evaluate its statistical characteristics (e.g. the chi-squared test for evaluating deviation from Poisson distribution). Due to the limited ability to directly count the number of encounters between a host cell and chytrid zoospores, we operationally regard the number of attached fungi as the number of successful infections.
Second, when we detect non-random interactions that deviate from a Poisson distribution, can we quantitatively describe the non-random infection patterns? Our solution is to develop a simple mathematical model that describes multiple encounter processes between host cells and chytrid zoospores. In the model, we assume that the encounter rate (and successful infection rate) depends on the number of infections that have already occurred. With a mathematical analysis to derive the equilibrium infection pattern, we will be able to derive a mathematical function for the infection number distributions, which includes both random and non-random interaction terms. Then, by fitting the mathematical function to the empirical patterns of multiple infections, we can statistically distinguish the infection patterns that emerged from the random encounter rate from those that could emerge from the encounter rate depending on the number of the existing infections.
Third, are there any chemical signals that can explain random and non-random interactions? Cell exudates are responsible for chemotaxis and host recognition by chytrid zoospores (Scholz et al. 2017). In terrestrial plants, herbivore in festation and disease infection induce compositional changes in volatile organic compounds (VOCs), which alters the behavior of enemies (predators, parasites, and disease agents) (Arimura et al. 2009), e.g. attracting towards or repelling from the infested/infected hosts. VOCs released from microbes (mVOCs) are also key drivers of below-ground interactions be tween plants and microbes and those between microbes (Schulz-Bohm et al. 2017). In aquatic systems, microalgae are also known to release VOCs and affect other organisms in the surrounding water (e.g. acting as the defense and foraging cues; Watson 2003, Fink 2007). However, no previous studies have examined whether and how chytrid infection changes VOC compositions from phytoplankton. Our strategy is to collect VOCs and compare their compositions between those from uninfected host cells and from infected cells. If there are VOCs specifically detected in uninfected or infected cells, these VOCs may act as (1) possible cues for chytrid zoospores to distinguish between un infected and infected host cells, and/or (2) possible proxies for host cell physiological status (e.g. resistance against or susceptibility to fungal infection). We hypothesize that such specificity of VOC composition exists for the host−chytrid pair that shows the non-random infection number distribution.
To answer these 3 questions, we compared multiple pairs of phytoplankton hosts and chytrid strains in infection and VOC experiments. Each pair consisted of a single host strain (from Ulnaria sp., Fragilaria crotonensis, or Staurastrum sp.) and a single chytrid strain. We then combined the statistical analysis on the experimental results and modeling analysis to characterize the interactions between phytoplankton hosts and chytrid parasites.
Cultures were grown in a temperature-and lightcontrolled incubator (IGB, Bender) at ~40 μE m −2 s −1 and 17°C. The phytoplankton strains (hereinafter, phyto-culture) and those infected by the corresponding chytrid strain (hereinafter, phyto−chytrid culture) were grown in 50 ml batch culture bottles with CHU-10 medium (Stein 1973). Every 2 wk, 10 ml from the phyto-culture bottle were inoculated with 2 ml from the phyto−chytrid culture bottle and 30 ml CHU-10 medium, in order to refresh the host−chytrid epidemic status. This is our regular protocol, which allows host cells to grow for 14 d until the chytrid causes the host population to crash (Senga et al. 2018). Therefore, we expected the newly produced chytrid zoospores to infect at least physiologically active host cells during our experiments, although we did not confirm this with growth curves. Both types of bottles (phytoculture and phyto−chytrid culture) were manually shaken every 1−2 d to reduce spatial heterogeneity in host and chytrid distribution within each bottle.

Experiments
After 4−18 d incubation (depending on the systems, see Table S1 in Supplement 2 at www.int-res. com/ articles/ suppl/a087p001_supp2.xlsx), we sampled 1 ml from the phyto−chytrid culture bottle and then fixed with 10 μl of Lugol's solution. For the Ulnaria systems, this experiment was done independently from the VOC experiment, whereas subsamples from the incubation for the VOC experiment (see also Section 2.4.1) were used for the Fragilaria and Staurastrum systems. Before counting the number of host cells and attached fungi, we added 25 μl of 3% Na 2 S 2 O 3 to remove the color from the Lugol's solution. Then, fungi attached to phytoplankton cells (attached zoospores and sporangia) were stained by adding 5 μl of 0.1% calcofluor white (Harrington & Hageage 2003). The number of host cells (in the Ulnaria and Staurastrum systems) or host colonies (Fragilaria system) with n attached fungi (n = 0, 1, 2, …) was counted at 400× magnification under a fluor escence microscope with both UV excitation (365 nm) and transmitted light. In order to exclude the effects of old infections from the inoculum of the phyto−chytrid culture bottle, we did not count the empty host cells (or colonies) and the empty sporangia. The prevalence of infection was calculated by dividing the number of infected cells (or colonies) (i.e. the cells with n > 0) with the total number of host cells (or colonies). In order to statistically distinguish Poisson and non-Poisson distributions, a substantial number of host cells should be counted. We counted at least 400 host cells (or colonies) for each phytoplankton−chytrid pair (Lund et al. 1958).

Overall strategy
Statistically, 3 types of infection number distribution are distinguishable. If zoospores do not discriminate uninfected and infected host cells, the en counter rate of zoospores on a host cell does not depend on the number of existing infections on the host cell, and it follows that the encounter between a host cell and zoospores occurs at random. In this case, the infection number distribution will follow a Poisson distribution (Anderson & May 1978). If the existing infection suppresses further infection through re pelling zoospores (or via the host's induced defense), the number of infections is likely to follow a uniform distribution. If the infection enhances further infection by attracting zoospores (or weakening the host's defense level), and/or if there are inherent variations between cells within a population in susceptibility to the infection, e.g. the variation in growth phase, it is likely to follow a clumped distribution (Anderson & May 1978). We used the chi-squared test to determine whether the infection number distribution follows a Poisson distribution or not. We also used the standardized Morisita index to check whether the distribution is random (following Poisson distribution), clumped, or uniform (Morisita 1962). In order to exclude the inherent heterogeneity in host cells and chytrid zoospores as much as possible, we use single-strain host and single-strain chytrid cultures, but under non-axenic conditions (e.g. bacteria were not excluded). We also prepared multiple pairs of host− chytrid cultures for comparison. given by, .

Chi-squared test
Using the function pchisq() in the package vegan (Oksanen et al. 2019), we calculated the p-value, i.e. the probability at which the chi-squared values become greater than the observed chi-squared value under the assumption of a Poisson distribution. When the p-value was small (p < 0.01), the null hypothesis, i.e. that multiple infections occur at random and the infection number follows a Poisson distribution, was rejected.

Morisita index
In order to quantify the clumpedness and uniformness, we used the standardized Morisita index of aggregation (Smith-Gill 1975), using the function dispindmorisita(). When the index value Imst < −0.5, the distribution is uniform, whereas the distribution is clumped when Imst > 0.5. When −0.5 ≤ Imst ≤ 0.5, the distribution is Poisson.

Model variables
We consider a host phytoplankton population that has cells in n states, whereby n represents the number of infections. The state n = 0 represents the uninfected state of cells, whereas n ≥ 1 represents the infected states of cells including multiple infections (n ≥ 2). Let H n (t) denote the abundance of cells with the state n at time t. We also set F(t) as the abundance of chytrid zoospores at time t. The following modeling follows the classical but standard framework for modeling aggregation in parasite− host systems (Adler & Kretzschmar 1992, Wilber et al. 2016a).

Local interactions between chytrid zoospores and host cells
We first defined the 'local' abundance of the existing number of infections on a single host cell as n. We then assumed the attachment of a single chytrid zoospore to host cells with n infections (its abundance is denoted as H n ), and the subsequent infection progress depends on n (Fig. 1). This new infection rate is given by a 0 exp(a 1 n+a 2 n 2 ). Here, the number of the existing infections (n) can be regarded as the local density of chytrid zoospores (or fungal cell structures at later infection stages) on a single host cell. The parameters a 0 and a 1 represent density-independent and density-dependent infection coefficients, respectively. The parameter a 2 is the coefficient for n 2 for representing a non-monotonic change of the new infection rate with the local density n when a 2 is unequal to zero. For example, with a 1 > 0 and a 2 < 0, the new infection rate initially increases with n but decrease when n is very large. During the progress of infection processes, the in fection may fail depending on the interactions between the newly infected zoospore and the existing matured infections. In this case, the host cell state shifts from n to n − 1, depending on the local infections (n), with the rate r 0 exp(r 1 n+r 2 n 2 ). Here, we operationally define it as the partial recovery rate where r 0 is the densityindependent recovery coefficient and r 1 and r 2 represent density-dependent re covery coefficients.
Temporal changes in H n with the infection and recovery processes (occur-ring at a daily scale; infectivity of zoospores is lost within a couple of days; Bruning 1991b) occur along with the demographic processes of phytoplankton hosts (e.g. doubling time falls into a weekly scale; Bruning 1991b) and new production and mortality of chytrid zoospores F. At the first step of developing a mathematical model for multiple infections in phytoplankton− chytrid interactions, we mathematically assumed that infection and recovery occur faster than the demographic processes (time-scale separation; Auger & Poggiale 1996) and we focused on the infection and recovery processes only. Then, for each cell state (n = 0, 1, 2, …), we can describe the dynamics of the state transition as: (1)

Infection number distribution
From the model assumption, for given total abundance of host cells ( ) and that of zoospores (F ) at time t, the frequency of different states of host cells (p n (t) = H n (t) / H T (t), n = 0, 1, 2, …), that is the infection number distribution, reaches equilibrium without the influence of demographic processes, as follows (see Text S1 in Supplement 1 at www.int-res. com/ articles/ suppl/a087p001_supp1.pdf): n n a n a n n a n a n n r n r n n r n r n n H H  (2) where α 2 and α 3 represent the net density-dependent effect, and is given as Also , and p 0 is a function of λ, α 2 , and α 3 . Note that this distribution is identical to the Poisson distribution, with p 0 = e −λ , when α 2 and α 3 are zero, whereas it deviates from the Poisson distribution when α 2 or α 3 is unequal to zero (also see Text S1 in Supplement 1).

Parameter estimation
We finally compared the goodness of fit (by Akaike information criterion, AIC) of the Poisson distribution and the theoretical distribution (Eq. 2) by fitting the theoretical distribution (Eq. 2) to the observed infection number distribution. With logarithmic conversion, we can apply a linear regression function in R (lm()) to the original nonlinear occurrence distribution function (from Eq. 2) using the following equation, where b subscript i (i = 0, 1, 2, 3) is the statistically fitted coefficient. We fit the Poisson distribution by using the estimated parameter , and 3 versions of the non-Poisson model (1 st order, 2 nd order, or 3 rd order forms with respect to n) (Eq. 3). Note that the Poisson fitting with a single parameter λ Ε and the linear form of Eq. (3) are potentially different because the latter has 2 parameters (b 0 , b 1 ). Only when the observations were significantly different from the Poisson expectation (tested by chi-squared test), did we compare the 3 versions of the non-Poisson models by the AIC() function in R.

Incubation
We used 3 phytoplankton hosts for our VOC experiments: Ulnaria sp., Staurastrum sp., and F. crotonensis. Among the multiple Ulnaria systems with different chytrid strains, we used the system with the C4 strain for the VOC experiments only. For the pre-incubation (starting at Day 0), we incubated 3 replicates of 200 ml bottles from the Ulnaria sp. and Staurastrum sp. phyto-culture, respectively. We also incubated 2 (or 1) bottles of 220 ml chytrid mixture, i.e. 200 ml Ulnaria sp. (or Staurastrum sp.) phyto-culture + 20 ml phyto− chytrid culture, for 5 d. For F. crotonensis, we incubated 6 replicates of 200 ml bottles from the Fragilaria sp. phyto-culture and added 30 ml of the chytrid mixture at Day 4 into 2 of these replicates, and thus we had 4 replicates of 8 d old host cultures in 200 ml bottles and 2 replicates of 5 d old chytrid mixtures in 230 ml bottles at the end of the pre-incubation (Day 8). The heterogeneity in pre-incubation conditions between the different host species originated from the fact that the pre-incubated samples were also used for the preliminary experiments to confirm whether it is possible to collect VOCs from this small culture volume (not shown). After the pre-incubation (at Day 5 for the Ulnaria and Staurastrum systems and at Day 8 for the Fragilaria system), we partly mixed the pre-incubated host culture, chytrid culture, and the fresh medium to prepare 200 ml bottles for hostonly treatments and host + chytrid treatment (Table S2 in Supplement 2). For every set-up, we prepared 3 replicates, and the incubation period until VOC collection was 9 d for the Ulnaria and Staurastrum systems (at Day 14) and 5 d for the Fragilaria system (at Day 13). We also continued the incubations of the leftover from the pre-incubations for counting the infection number for both the Staurastrum (at Day 18) and Fragilaria systems (at Day 13) (Table S1).

VOC collection
After incubation, the host cell density was estimated by chlorophyll a fluorescence and the calibration curves between counted cell density and fluorescence, which was measured with a fluorometer (Laboratory Fluorometer with chlorophyll a in vivo module, Trilogy ® ). The estimated cell density (×10 3 cells ml −1 ) was 36 ± 20, 5.5 ± 3.3, and 89 ± 20 for the Ulnaria, Staurastrum, and Fragilaria systems, respectively (average ± SD). Then, 200 ml of each culture were concentrated into 15 ml by using a 10 μm nylon mesh. The concentrated culture was then transferred to a 30 ml glass vial and 1 Twister ® (GERSTEL) was placed in the culture and 2 more twisters were set in the headspace of the culture (Fig. 2) during the collection of volatiles. Twisters were allowed to absorb the volatiles dissolved in the culture media and those released into the headspace for 24 h. In order to distinguish VOCs that were naturally dissolved in the media and those from the cultured organisms, a single control vial with 15 ml WC medium was prepared with the same twister placement.

VOC analysis
The collected volatiles were directly desorbed in the injector with a thermal desorption system, a cooled injection system, and cold trap system (Gerstel) and then analyzed by using a GC-MS (6890/5972 GC-MS system, (Agilent Technologies) with an HP-5 MS capillary column (length, 30 m; inner diam., 0.25 mm; film thickness, 0.25 μm). The oven temperature of the GC was programmed to rise from 40°C (9 min hold) to 280°C at 5°C min −1 . Headspace volatiles were identified by comparing their mass spectra to those in the Wiley7N database.

Statistical analysis
In order to focus on VOCs that were released from the organisms, which included phytoplankton host and chytrid as well as bacteria, the shared components between the treatments and control were excluded for further statistical analyses, except for limonene. Although limonene was also detected in the control, the difference in magnitude between the control and treatments were substantial so that the values from each treatment minus the value from the control were used for the analysis.
We used R version 3.6.2 and the vegan package for the following statistical analysis. We used Bray-Curtis dissimilarity to evaluate differences in the overall VOC composition between treatments and replicates. For the visualization of the difference, we used a principal coordinate analysis (PCoA) plot by the function cmdscale(). To statistically test the difference in composition, we used permutation MANOVA by the functions adonis() and capscale(), checking the heterogeneity in variance via PERMDISP using the function betadisper(). We also identified major VOCs that are responsible for the differences between treatments via SIMPER analysis using the function simper(). For graphics, we used the ggplot2 package. For multiple comparisons, we used Bonferroni adjustment.

Infection number distribution and its statistical characteristics
We found that the infections in the Ulnaria system with the C4 strain and those in the Fragilaria system deviated from a Poisson distribution (chi-squared test;  Table 1). Numbers of the counted host cells and infections, as well as the prevalence of infections, were comparable between the 3 systems (Table 1). However, host cell density was much lower in the Staurastrum system than the other 2 systems. In fact, other Ulnaria systems also showed a clumped distribution (Table S3 in Supplement 2).

Fitting to the mathematical model
By comparing the AIC values, we found that the 2 nd order model (b 2 ≠ 0, b 3 = 0 in Eq. 3) and the 3 rd order model (b 3 ≠ 0) is the best model for the Ulnaria system (with the C4 strain chytrid) and the Fragilaria system, respectively (Fig. 3a,b). With the other strains (C3, C10, D11, and F12), Ulnaria's infection number distribution was best explained by either the 2 nd or 3 rd order model (Fig. S1). The fitting of the proposed model equations (Eq. 3) into the result from the Staurastrum system ( Fig. 3c) is just a reference, since its distribution did not significantly deviate from a Poisson distribution (Table 1). The estimate of parameters, adjusted R 2 values, and AIC values are summarized in Table S4 in Supplement 2, and the results from the other Ulnaria systems are shown in Fig. S1.

VOC composition
In total, 13 VOCs were detected from uninfected phytoplankton and infected phytoplankton in total, which were not detected in the control (growth media only), or their amount in the control was much lower (only for limonene) ( Table S5 in Supplement 2). The compositional variations in VOCs were visualized by a PCoA plot, and the first 2 axes explained 33.7% and 23.4% of the whole variations (Fig. 4). The contribution of each VOC component to these first 2 axes is summarized in Table S6 in Supplement 2. Overall VOC composition was statistically different between the Ulnaria and the Fragilaria systems (pairwise PERMANOVA, constrained variation = 0.42, p = 0.003) and between the Ulnaria and the Staurastrum systems (pairwise PERMANOVA, constrained variation = 0.50, p = 0.0005), but the difference between the Fragilaria and the Staurastrum systems was not significant (pairwise PERMANOVA, constrained variation = 0.24, p = 0.11) (Fig. 4). Since dispersion was not different between species (PER -DISP, p = 0.377), the difference detected in PERM-ANOVA represents the averaged compositional differences. On the other hand, the overall VOC composition was not different between uninfected and in fected treatments (PERMANOVA, p = 0.365).
More specifically, SIMPER analysis found 6 representative compounds to distinguish species (limonene, octyl ether, phytane, unidentified component, 6-dode- 3), respectively; blue and light-blue lines: the 2 nd order (b 2 > 0, b 3 = 0) and 3 rd order (b 3 > 0) model fitting, respectively. Note the different x-axis scales cenol, and beta-ionone), which, except for 6-dodecenol, were also the representatives to distinguish chytrid infection (Fig. 5). In particular, the levels of beta-ionone in the uninfected Staurastrum system were significantly higher than in other species, but chytrid infection significantly reduced it to levels seen in the other cultures (linear model, p = 0.00173). The unidentified VOC, which is close to (E)-10-Heptadecen-8-ynoic acid methyl ester, and dodecenol were detected from the Ulnaria system only, whereas phytane was not detected from the Staurastrum system. Phytane apparently increased with chytrid infection in the Fragilaria system, but the difference was not statistically significant.

Occurrence of individual-level heterogeneity and its non-randomness
With a simple microscopic observation to obtain the distribution of infection number, we successfully demonstrate the occurrence of individual-level heterogeneity in terms of infection number and that it is possible to distinguish random and non-random het-erogeneity. A clumped chytrid parasite distribution on Ulnaria sp. and Fragilaria crotonensis implies 2 potential mechanisms to realize the non-randomness, which are not mutually exclusive. First, it could result from the presence of internal physiological heterogeneity between uninfected host individuals in the population (see Section 4.4). Alternatively, chytrid infection further enhanced the subsequent re-infection rate, which resulted in a higher frequency of host cells with many at tached fungi than expected from a Poisson distribution. This possibility was further examined by a mathematical model (see the next section).

Modeling of transition of host cell states regarding the number of infections
The mathematical model predicted a weighted Poisson distribution (WPD; Eq. 2). To represent over-and underdispersion, various WPD have been proposed, such as a 2-parameter exponentially weighted Poisson distribution (EWP) (Sellers et al. 2012). In the present study, we derived another EWP with the novel weight function exp(α 2 n 2 +α 3 n 3 ). When a chi-squared test suggested that the fitting of the Poisson distribution was not good, we found that nonlinear models (in the 2 nd order or 3 rd order) fitted better than the linear model (1 st order) that assumes constant infection and partial recovery rates. This result indicates that infection and/or partial recovery could occur, depending on the number of existing infections on each host cell.

Characteristics of VOCs and their linkages to infection number distribution
Based on the implications from the non-random infection number distribution, we hypothesized that there should be a difference in VOC composition between uninfected and infected host populations for the Ulnaria and Fragilaria systems, but not for the Staurastrum system. This is based on the assumption that the altered infection rate by the existing infections resulted from the infection-specific VOC composition, acting as infochemicals (either attractants or   cally no differences in the overall VOC composition between uninfected and infected populations for all 3 tested host species. At the same time, the second hypothesis was supported by the lack of any VOCs with antifungal effects enhanced by chytrid infections. However, there was a statistically significant difference in VOC composition between the host species, and a deeper focus on the 6 representative compounds gave us some insights into the ecological roles of VOCs in phytoplankton−chytrid interactions. First, we identified 3 species-specific compounds: 6dodecenol and an unidentified one (its signal was close to (E)-10-Heptadecen-8-ynoic acid methyl ester, but the quality was low) were released only from Ulna ria sp., and phytane was not released from Staurastrum sp. 6-Dodecenol has been detected from the bacterium Pseudomonas fragi (Ercolini et al. 2009). (E)-10-Heptadecen-8-ynoic acid methyl ester has been detected from the green alga Dunaliella salina (Stocker et al. 2008, Atasever-Arslan et al. 2015. Phytane (2, 6,10,14-tetramethyl hexadecane) is formed from phytol, a constituent of chlorophyll, which has been detected from bacteria (Rontani & Volkman 2003). Although the ecological functions of these 3 compounds remain unclear, they could act as infochemicals released from host cells and/or attached bacteria, allowing chytrid zoospores to distinguish between host and non-host phytoplankton species.
Second, limonene, a monoterpene, which is produced by many plants and algae including diatoms (Wise 2003, Meskhidze et al. 2015 and cyanobacteria, has allelopathic effects on other phytoplankton (Hu et al. 2014, Zuo 2019, and has antifungal (Ünal et al. 2012) andantibacterial (Pathirana et al. 2018) effects. Although the amount of limonene was not significantly different between host species and between infection status, it could be involved in the constitutive defense of hosts to fungal/bacterial infections. Beta-ionone, a sesquiterpene, which is produced by many species of bacteria (from the mVOC database: Lemfack et al. 2018), fungi (Ziegenbein et al. 2006), cyanobacteria (Jüttner 1984), and green algae (Lafarge & Cayot 2019), exerts negative effects on green algal growth (Ikawa et al. 2001) and potentially on fungi (Abdel-Hafez et al. 2015). Interestingly, the amount of the release of beta-ionone was the largest from uninfected Staurastrum sp. Despite a lower population density, the amount of this VOC was reduced after chytrid infection (Fig. 5). Since the prevalence of chytrid infection on Staurastrum sp. was not lower than on the 2 other host species (Table 1), it is relatively unlikely that beta-ionone is involved in resistance against chytrid infection.
Third, octyl ester (di-octyl ether) has been reported as a VOC released from wood-decaying fungi in soil (O'Leary et al. 2019). Therefore, it could originate from the chytrid fungi in our systems. However, octyl ester was also detected from the uninfected host populations. Since our culture systems were non-axenic, the source of octyl ester could be also bacteria in the surrounding water. We cannot rule out the possibility that other fungi in the bottles released octyl ester, yet we never observed growth of other fungi in the cultures. Also, it could stem from the host itself, since it is also reported that terrestrial plants, liverworts, release this compound (Ludwiczuk et al. 2008). Due to the data limitation, it is not clear if it was released from the hosts, chytrid fungi, or other fungi and bacteria in the surrounding water. It is reported that Acinetobacter can metabolize it (Modrzakowski et al. 1977), yet, its effects on phytoplankton and fungi remain fully unexplored. It will be necessary to improve the filtering method after incubations for clearly separating VOCs from phytoplankton− chytrids and those from bacteria.

Potential mechanisms that differentiate random and non-random distributions
For the Ulnaria and Fragilaria systems, the overdispersion suggested by the chi-squared test and a good fit to the negative binomial fitting (not shown) indicate that such distribution can emerge as the superposition of many Poisson distributions with different λ (=averaged infection number) values. It biologically implies that the host population consisted of many subpopulations with different physiological status or at different cell cycle stages. As a result, susceptibility to zoospore attraction or subsequent infection processes differed between them. Such variations include the asynchrony in cell cycle stages, cell age differences partly originating from our incubation settings, and unknown inherent phenotypic variations (e.g. growth rate) maintained by epigenetics (Traller & Hildebrand 2013, Tirichine et al. 2017. If the latter is the case, it implies the presence of subgroups with a high susceptibility to chytrid infection even within a genetically uniform population. Such phenotypic variations might have emerged due to genetic heterogeneity, since phytoplankton strains were isolated a couple of years before the experiments and some mutations might have accumulated. From a morphological and physiological point of view, however, individual cells behaved similarly to each other. Alternatively, but not mutually exclusively, a clumped distribution could emerge from densitydependent sequential encounter and attachment of chytrid zoospores to the host cells. Since we did not find any evidence for VOCs being actively involved in differentiating uninfected and infected cells, it is likely that zoospores rely on other types of chemical cues released from host cells. Other studies have suggested that active release or passive leakage of exudates during photosynthesis provide chemical cues for zoospores to find their host (Scholz et al. 2017). The infection may also change the host's cell physiology and the activity of attached bacteria, resulting in a change in the dissolved organic matter (DOM) compositions released from the host cells (Senga et al. 2018), which may act as potential chemical cues for zoospores.
The Poisson distribution in our Staurastrum system indicates that there were neither uninfected subpopulations with different susceptibility to parasitism nor differences in susceptibility between uninfected and infected cells. This finding is in contrast to the commonly observed clumped distribution of parasites in other systems (Crofton 1971, Poulin 2007. Theoretically speaking, the infection number distribution of phytoplankton is likely to deviate from random, falling into the clumped distribution confirmed in both the Ulnaria and Fragilaria systems, or a uniform distribution when the induced resistance of the host cells and/or resource competition between attached fungi are strong (Anderson & May 1978, May & Anderson 1978, Müller-Graf et al. 2001). In addition, the absence of a non-random parasite distribution indicates that interactions between attached fungi on each host cell, which enhances or reduces the infection rate of infected cells (leading to clumped or uniform distribution, respectively), were negligible in the Staurastrum system. Although it would be possible that Staurastrum sp. has some inducible mechanisms against any fungal infections (cf. induced resistance against zooplankton grazing; Wiltshire et al. 2003), its lower average number of infections (λ E = 0.497) and the lower maximum number of infections per cell (=4) than in the other 2 systems could explain the lack of substantial interactions between attached fungi.
There are also other potential mechanisms for the variation in the infection number distribution (random vs. clumped) and VOC compositions, which arise from differences in biology and incubation conditions between species.
First, differences in the reproductive mode of chytrid fungi may affect the degree of chytrid aggre-gation. In sexually reproducing organisms, intermediate levels of aggregation may be favored by natural selection as it increases the probability of mating (Shaw & Dobson 1995). Chytrids are known to produce resting spores which allow their survival in times of host absence or adverse environmental conditions. Resting spore formation occurs via sexual reproduction on the host (e.g. conjugation of 2 chytrid thalli on the same host cell) (Doggett & Porter 1996). In this case a certain level of infection aggregation is necessary to assure mating success and thus resting spore formation. However, there are also chytrid species that lack sexual reproduction and can produce resting spores asexually, as is the case for the chytrid species infecting Staurastrum in this study (Van den Wyngaert et al. 2017). Such an asexual reproductive strategy, which has no need for nearby conspecifics, may cause less selection pressure towards aggregation, and potentially explains the random distribution observed for this chytrid.
Second, differences in the effects of chytrids on host population dynamics may affect the chytrid distribution. Although we assumed in the model that the distribution pattern and population dynamics (replication and mortality) occur at different timescales, population dynamics and demographic processes (e.g. infection-induced host mortality) are also im portant (Anderson & Gordon 1982, McCallum & Ander son 1984, Wilber et al. 2016b. For example, if the lower maximum infection number (=4) in the Stau rastrum system compared to the other 2 species was due to higher infection-induced host mortality (when n ≥ 5) and/or lower quality of hosts that prevented further multiple infections, it may contribute to lower dispersion of this system than clumped distribution (Anderson & Gordon 1982).
The differences in the best model selected and estimated parameter values (Tables S3 & S4) could result from both chytrid strain-specific effects and differences in incubation time (and thus difference in the prevalence of infection). Nevertheless, the results from the Ulnaria system clearly indicate that the clumped distribution of chytrids on Ulnaria sp. is a robust pattern and is independent of the used chytrid strains and incubation time. For the Staurastrum system, a lower cell density and the lower growth rate due to the longer incubation period (although we did not directly investigate the growth curve) than for the other 2 host species systems may have had a negative effect on the encounter rate with chytrids. On the other hand, its greater cell size (biovolume) could have had a positive effect on the encounter rate. Therefore, the net effects of these contrasting char-acteristics on the number of multiple infections and its random infection distribution pattern remain unclear. Similarly, differences in VOC composition be tween host species and between uninfected and infected cells could partly result from the differences in the incubation period (Table S2).

Potential impact of clumped distribution on population dynamics
Non-randomness of the infection number distribution and the degree of the aggregation will affect the reproduction number (R), which is a core index to evaluate the potential of any disease spread. This is defined as the total number of new infections caused by a single infective agent (May 2006, Sofonea et al. 2015. If the burst size of the chytrid (=the number of newly produced zoospores per host cell) linearly increases with the infection number n, the multiple infections will substantially contribute to R. However, since the resource for chytrid growth and reproduction within a host cell is finite, it is likely that the reproductive outcome per infection decreases with the infection number n due to resource partitioning between multiple attached chytrids (i.e. within-host competition ;Bruning 1991b). If this is the case, the relationship between the burst size and n is a critical determinant of R. In order to better elucidate this point, it will be necessary to directly evaluate the individual-level heterogeneity in newly produced zoospores and its dependence on the degree of multiple infections in a future study.

Concluding remarks and outlook
Our study has demonstrated the effectiveness of mathematical modeling and chemical ecology approaches to investigate individual-level heterogeneity in phytoplankton−chytrid interactions. These approaches indicate multiple potential mechanisms that drive individual-level interactions between host cells and chytrids, but which are not mutually exclusive, and the exact mechanisms still remain unclear. In order to directly compare these 2 mechanisms -(1) internal physiological heterogeneity between uninfected host individuals, and (2) enhancement of subsequent re-infection rate -we need to track the se quence of attachment of multiple zoospores on a single host cell. With such sequential data over time, we will be able to directly quantify the dependence of the timing of infections and recovery on the num-ber of parasite infections, without assuming equilibrium of host−chytrid infection processes. The characterization of zoospore swimming patterns (e.g. random or directed) and their response to species-specific VOCs will also be useful to make clearer the interspecific variations in infection distribution patterns. Direct observation and quantification of individual heterogeneity in other epidemiological parameters (e.g. burst size and development time) using culture systems as well as quantifying the heterogeneity in natural systems will also be necessary to narrow down the proposed mechanisms and hypotheses. In addition, further development of a mathematical model that couples multiple infection processes and population dynamics (Adler & Kretzschmar 1992, Cosner et al. 1999, Wilber et al. 2016a, and further exploring chemical ecology in chytrid parasitism, will be the next steps for a better understanding of chytrid infection dynamics.