Droplet digital PCR as a tool for investigating dynamics of cryptic symbionts

Abstract Interactions among symbiotic organisms and their hosts are major drivers of ecological and evolutionary processes. Monitoring the infection patterns among natural populations and identifying factors affecting these interactions are critical for understanding symbiont–host relationships. However, many of these interactions remain understudied since the knowledge about the symbiont species is lacking, which hinders the development of appropriate tools. In this study, we developed a digital droplet PCR (ddPCR) assay based on apicomplexan COX1 gene to detect an undescribed agamococcidian symbiont. We show that the method gives precise and reproducible results and enables detecting cryptic symbionts in low target concentration. We further exemplify the assay's use to survey seasonally sampled natural host (Pygospio elegans) populations for symbiont infection dynamics. We found that symbiont prevalence differs spatially but does not show seasonal changes. Infection load differed between populations and was low in spring and significantly increased towards fall in all populations. We also found that the symbiont prevalence is affected by host length and population density. Larger hosts were more likely to be infected, and high host densities were found to have a lower probability of infection. The observed variations could be due to characteristics of both symbiont and host biology, especially the seasonal variation in encounter rates. Our findings show that the developed ddPCR assay is a robust tool for detecting undescribed symbionts that are otherwise difficult to quantify, enabling further insight into the impact cryptic symbionts have on their hosts.


| INTRODUC TI ON
Interactions between symbiotic organisms and their hosts dramatically influence not only organismal ecology and evolution, but also the dynamics of entire ecosystems (Faust & Raes, 2012;Godfrey-Smith, 2015).
For the majority of symbioses, however, we do not know how the interacting species affect each other, much less the broader consequences of the symbiosis for communities or ecosystems. Often, it is the lack of tools for effectively monitoring or investigating symbiotic interactions that prohibits progress. In addition, a lack of general knowledge about the symbiotic organisms can limit the development of these tools (Pawlowski et al., 2012), leaving much of the functional and specieslevel diversity understudied.
Apicomplexans are a diverse group of microbial eukaryotes that include some of the most studied parasites with significant social and economic consequences (e.g., the genera Plasmodium, Toxoplasma, Eimeria, and Cryptosporidium) (Seeber & Steinfelder, 2016). Infecting majority of vertebrates and invertebrates in terrestrial and aquatic environments, apicomplexan species diversity is estimated to be over a million (Pawlowski et al., 2012). However, the bulk of apicomplexan diversity, especially in marine environments, is largely undescribed and contains many cryptic species (Janouškovec et al., 2015;Xavier et al., 2018). Additionally, despite the fact that some groups are known to be common and prevalent symbionts of marine invertebrates, little is known about the nature of their interactions (mutualistic or parasitic) with their hosts (Rueckert et al., 2019). Studying symbiont dynamics, even when the symbiont is cryptic, could help to shed light on the nature of their interaction with the host and their effect on host populations.
The presence of apicomplexans within host individuals and populations is traditionally determined with microscopy, but this can be difficult and labor intensive when studying large samples and may overlook small or cryptic microbial eukaryotes. In such cases, the use of molecular tools, such as PCR methods or amplicon sequencing, often are more appropriate. For instance, real-time quantitative PCR (qPCR) has been widely used in detecting and quantifying protistan parasites (e.g., Maia et al., 2014). Digital PCR (dPCR) is a third generation PCR method for quantification of target molecules in a nucleic acid sample (Li et al., 2018). It is based on partitioning and randomly distributing the sample into small partitions before PCR amplification, which takes place in each partition separately (Hindson et al., 2011). After amplification, the end-point reaction is visualized to determine the fraction of positive partitions (Hindson et al., 2011(Hindson et al., , 2013. From this fraction, the concentration of the target can be estimated using Poisson statistics (Sykeset al., 1992). In droplet digital PCR (ddPCR), the sample is partitioned into up to 20,000-nl-sized droplets by water-oil emulsion and the target concentration is therefore determined from the fraction of positive droplets. The advantage of partitioning lies in increased resolution and sensitivity; hence, ddPCR applications have been widely used in clinical studies, GMO detection, and food security (Demeke & Dobnik, 2018;McMahon et al., 2017;Sedlak et al., 2014).
Compared with traditional real-time qPCR, ddPCR has several major advantages. For example, it has higher precision and dayto-day reproducibility, especially when targeting rare copies (Brys et al., 2020;Doi et al., 2015;Hindson et al., 2013). ddPCR assays also tolerate inhibition better than qPCR (Dingle et al., 2013;Poh et al., 2020;Rački et al., 2014), and the end-point measurement of the nucleic acid quantitation is not affected by amplification efficiency. In addition, ddPCR does not require standard curves (Volgenstein & Kinzler, 1999). Therefore, the concentration of the target copies is an absolute measurement that is not reliant on Cq values. Furthermore, there is increased statistical power of ddPCR over qPCR (Taylor et al., 2017). These advantages make the method particularly applicable for molecular identification of symbionts. The increase of studies using ddPCR for quantifying haemoprotozoan infections (e.g., Plasmodium and Babesia) (Koepfli et al., 2016;Srisutham et al., 2017;Wilsonet al., 2015) and for detecting parasites from environmental samples (Mulero et al., 2020;Rusch et al., 2018) has already shown that the method is robust, reproducible, and capable of detecting rare copies of target DNA.
Despite the aforementioned advantages, ddPCR has not been used to its full potential, in particular for investigating cryptic symbionts and their ecological interactions with hosts. In this study, we show that ddPCR can shed light on biological interactions with precise estimates of cryptic, unculturable symbionts in natural host populations, indicating that the method is useful in molecular ecology research. We introduce a detailed ddPCR protocol and demonstrate its use by investigating the seasonal dynamics of an undescribed cryptic agamococcidian (Apicomplexa) infection for the first time in four polychaete (Pygospio elegans; Claparède, 1863) populations. By using specific primers, our assay targets the mitochondrial cytochrome oxidase c subunit 1 gene (COX1) of the agamococcidian symbiont and can be used effectively to determine symbiont prevalence and infection load from whole host DNA extracts.

| Study organisms
The host species, P. elegans, is a small polychaete worm that inhabits sandy coastal habitats throughout the northern hemisphere.
This species can be a dominant member of benthic communities (Bolam, 2004;Bolam & Fernandes, 2003) and is an important prey item for other invertebrates and fish (Mattila, 1997). Pygospio elegans is known to host at least two apicomplexan symbionts: the archigregarine Selenidium pygospionis (Paskerova et al., 2018) and the eugregarine Polyrhabdina pygospionis (Paskerova et al., 2021), which both inhabit the worm's intestine. In our previous study of the host's transcriptome (Heikkinen et al., 2017), we detected the presence of a third apicomplexan symbiont. Based on 18S rDNA sequence similarity, this symbiont is most likely an agamococcidian (Order Agamococcidiorida; Levine, 1976), but its location within the host and definitive identification is not yet known. Agamococcidians are a small group of coccidian-like symbionts currently comprising two monogeneric families Gemmocystidae and Rhytidocystidae (Levine, 1976;Upton & Peters, 1986), but the taxonomy of the group has been questioned (Janouškovec et al., 2019;Mathur et al., 2020).
Rhytidocystidae includes species infecting midguts of marine polychaetes (Leander & Ramey, 2006;Miroliubova et al., 2020;Rueckert & Leander, 2009), and some have been suggested to have intracellular life stages at least at an early stage of development within the hosts (Miroliubova et al., 2020). Because research on this group has mostly focused on resolving their phylogenetic position and describing the species, the nature of their interaction with their hosts has not been studied previously.

| ddPCR assay specificity
Our goal was to produce an assay that could be used to detect and estimate infection loads of a cryptic symbiont from total DNA extracts (containing a mixture of host and symbiont DNA). A limited dataset of potential genes was available from our previous transcriptome study of the host (Heikkinen et al., 2017). We chose to design primers to amplify the presumed mitochondrial COX1 gene of the undescribed agamococcidian, as it showed sufficient divergence from the host's COX1 sequence. Primer3 software (https:// prime r3.ut.ee) was used to choose appropriate primers: forward primer ApiCox1F (5′-ACT GGT CTA TCA AGT GTA CTG GC-3′) and reverse primer ApiCox1R (5′-GAT CAC CAC TAA ATT CAG GGT CA-3′) to amplify 226 bp of the target gene. Amplicons were identical in sequence to the previously obtained transcript and showed high similarity to transcripts obtained from Rhytidocystis sp. ex.
Non-attached gregarines were hand-picked from dissected host intestines by micromanipulation and washed three times with filtered (Millipore 0.2 µm) seawater. The cells were then transferred into microtubes, pelleted down with centrifugation to remove excess water, and stored in 100% EtOH for transport to our laboratory. Prior to DNA extraction, EtOH was evaporated from the samples and the DNA was extracted using MasterPure™ Complete DNA and RNA Purification Kit (Lucigen) following the manufacturer's protocol, except for the elution step, which was done in 10 μl of EB buffer. DNA concentration was checked using a Qubit 4.0 Fluorometer with 1X dsDNA HS Assay (Thermo Fisher Scientific).
2.3 | ddPCR protocol: reproducibility, dynamic range, and the limit of detection ddPCR was performed using Bio-Rad's QX200™ Droplet Digital™ PCR System. Reproducibility of the assay was inspected with four replicates of three DNA samples: host with high (approx. 100 copies/µl), intermediate (approx. 50 copies/µl), and low (approx. 2 copies/µl) infection load. Because samples are partitioned into droplets in sets of eight in the QX200 system, two replicates of each sample were used in two different droplet generation events (four replicates for each sample). To inspect the dynamic range of the assay and the limit of detection (LOD), a 10-fold dilution series (1, 1:10, 1:100, and 1:1000) of the DNA extracted from the host with high infection load was performed in 10 replicates. Linearity over the dynamic range was determined by the coefficient of correlation r 2 , calculated on the target concentration (copies/µl) measured in the dilution series replicates. Assay repeatability was determined by the % coefficient of variation (%CV = concentration standard deviation/concentration mean * 100) between the replicates. LOD is defined as the lowest target copy number in a sample that can reliably be detected. In this study, LOD was determined as the lowest concentration level for which all 10 replicates resulted in at least three positive droplets per replicate. In all experiments, a negative control containing nucleotide free water was used. The reaction mix was prepared to a volume of 20 µl per sample. We used 2X QX200™ ddPCR™ EvaGreen ® (Bio-Rad) reagent mix. Primers were added to the mix in 1 µM together with 4.6 µl of sterile water. The reaction mix was divided into individual 0.5 ml microfuge tubes, and 2 µl of the DNA templates (varying concentration) was added to each tube so that the final reaction volume was 22 µl. Samples were partitioned into nanoliter-sized droplets with the QX200 Droplet Generator (Bio-Rad) using singleuse DG8 cartridges and Droplet Generation Oil (Bio-Rad). Twenty microliters of each reaction mix was loaded to the cartridge, and the emulsion was made with 70 µl of oil. The resulting droplets were manually transferred with a multichannel pipet to a ddPCR™ 96-well PCR plate (Bio-Rad), which was heat-sealed with a foil cover.
The droplets were then subjected to thermocycling using a Bio-Rad C1000 thermocycler with a ramp rate of 2°C/s in each step.
Initial denaturation of the DNA was done at 95°C for 3 min, after which the denaturation, primer annealing, and target extension steps were repeated for 40 cycles. The denaturation step was done at 95°C for 30 s, annealing temperature for the primers was optimized at 58°C for 1 min, and the target extension step was done at 72°C for 2 min. After the cycles, a signal stabilization step from 5 min at 4°C to 5 min at 90°C was added. Following the amplification, the droplets were immediately read with Bio-Rad's Droplet Reader.

| ddPCR data analysis
Absolute quantification of target gene copies was done with default ABS settings in QuantaSoft Analysis Pro 2.0 software (Bio-Rad). The ABS experiment estimates the concentration of the target in copies per microliter of the final 1X ddPCR reaction. Because one droplet can harbor one or more copies of the target, or none, the target concentration (copies/µl) is estimated by the software by calculating the mean copies per partition (λ) following Equation (1), where n is the total number of accepted droplets and k is the number of positive droplets counted.
(1) = − ln 1 − k n The infection load (copies/ng total DNA) was calculated using the reaction mix volume (22 µl), the sample volume (2 µl), and the concentration (ng/µl) of the DNA sample following Equation (2), where C ng is the number of copies per nanogram of total DNA, C ddPCR is the reaction concentration (copies/µl) given by QuantaSoft Analysis Pro, V r is the reaction mix volume, V s is the sample volume, and C DNA is the concentration (ng/µl) of total DNA sample.
Only reactions that had ≥10,000 droplets were included in further analyses. A threshold to separate the target positive and negative droplets was manually set in relation to the negative control by visual inspection.

| Survey of host populations
To to February 2015 to describe the reproductive dynamics of the host (see Thonig et al., 2016, for sampling details). In addition, the authors determined size (length in micrometer from eyespot to the beginning of gills) for each host and population density (individuals in square meter). Environmental parameters (salinity, temperature, and sediment characteristics) were also collected seasonally (see Thonig et al., 2016).

DNA was extracted from whole host individuals using DNeasy
Blood and Tissue Kit (Qiagen) following the manufacturer's protocol for animal tissue and stored in −20°C (Thonig et al., 2017), and DNA concentration was measured with Qubit 4.0 Fluorometer with 1X dsDNA HS Assay (Thermo Fisher Scientific). In the current study, a subset of those samples was analyzed with ddPCR (Table 1). These sampling times were chosen for quantification of agamococcidian infection due to population specific changes in host reproduction and density patterns observed in these months (Thonig et al., 2016).
We used the ddPCR protocol described earlier to detect the cryptic symbiont (prevalence) and estimate its abundance (infection load). Prevalence of the symbiont in the host population was measured as the proportion of infected P. elegans individuals in the sampled population. An individual host was considered infected if more than 0 copies of the target gene per ng of total DNA were detected with the ddPCR assay. Infection load was defined as the number of COX1 gene copies/ng of total DNA, as in Equation (2), and analyzed with only infected hosts. Aggregation of the infection was inspected by calculating variance to mean ratios. (2) Location of the four sampling sites in Isefjord-Roskilde Fjord complex, Denmark

| Statistical analyses
In the ddPCR assay reproducibility experiment, the normality of the data was checked with Shapiro-Wilk test. The concentration did not follow a normal distribution; hence, a non-parametric Wilcoxon rank sum test was used to compare the replicates of high, intermediate, and low infection load medians between the two droplet generation events. The linearity of dynamic range was determined using least squares regression analysis of the linear relationship between sample quantity (fold) and COX1 concentration (copies/µl) in each dilution series replicates. The correlation of droplet count and the absolute concentration was analyzed for all replicates used in the two experiments described earlier (n = 51) with Pearson's correlation.
To avoid multicollinearity issues caused by significant correla-  iterations per imputation. The parameters were estimated in each imputed dataset separately and combined using Rubin's rules (Rubin, 1987   no significant differences in the COX1 concentration between the droplet generation events were detected (Wilcoxon rank sum test: W = 18, p = 1.000), the median shift being 0.27 copies/µl (Figure 2).
Regression analysis showed that the concentration of the 10fold dilution series had good linearity (R 2 > .99, p < .001) (Figure 4).
However, one replicate did not show any amplification; hence, LOD was determined to be 0.86 COX1 copies/µl (1:100 dilution), as all 10 replicates had more than 3 positive droplets (Table 4). %CV was lowest for the undiluted original DNA sample and highest in 1:1000 dilution (Table 4).
Absolute concentration measurements were independent of the droplet count (Pearson's correlation: t = 0.449, df =49, p = .655) ( Figure 3). The droplet count varied from 10,239 to 17,553. Positive and negative droplets were easily separated from each other with threshold set to amplitude 13,033, mean amplitude for positives being 16,788, and negatives 5028. Some "rain" between the positive and negative clusters was observed, potentially because of the formation of primer dimer and the presence of background amplification (Figures 2a and 4a). The ddPCR assay did not amplify other apicomplexans (S. pygospionis and P. pygospionis) possibly present in the total DNA extractions from hosts (concentration 0 copies/µl) ( Figure 2a).

| Prevalence of infection within host populations
Prevalence of infection differed between the populations (Χ 2 = 64.9, df =3, p > .001). The highest proportion of infected hosts ranged from 85% to 100% in Herslev, where it was 94.9% more likely that the worms were infected than in Lammefjord, 77.1% more likely than in Vellerup, and 91.8% than in Lynaes (Table 5). The lowest prevalence was found in Lammefjord, where fewer than 50% of hosts were infected throughout the sampling period ( Figure 4, Table 5). No clear seasonal pattern common to the four populations was found (Χ 2 = 4.3, df =4, p = .36). The highest prevalence was observed in F I G U R E 3 Non-significant Pearson's correlation between the COX1 concentration (copies/µl) and the droplet count (n = 51) F I G U R E 4 Ten-fold dilution series to detect the lowest limit for agamococcidian detection for the ddPCR assay. (a) 1D ddPCR plot. x-axis shows an example of individual reactions for the undiluted, 1:10, 1:100, and 1:1000 diluted replicates, and y-axis shows the amplitude of the fluorescence. Positive reactions (blue) and negative reactions (grey) are separated by manually set threshold (pink line, at amplitude 11,309). Negative controls (neg) did not show any amplification above the threshold. (b) Dilution series sample quantity in log scale plotted against COX1 concentration (copies/μl) in log scale in each replicate. Significant correlation coefficient (r 2 = .99, p < .001) between the absolute concentration (COX1 copies/µl) and 10-fold dilution series shows good linearity of the assay Herslev in August (of all studied populations and sampling times) and the lowest in Vellerup in May ( Figure 5). Although the probability of being infected increased towards November, the increase was not significant (Table 5).
Larger worms had an increased probability to be infected, but the effect was small: For every micrometer increase in length, the probability of being infected was 0.1% higher (Table 5). When host population densities were high, the prevalence of infection was significantly lower, but the effect was small (Table 5). For every host individual increase within a square meter, the probability of infection decreased to 0.04%. Organic content (%) in the sediment was found to decrease the proportion of infected hosts by 67%, but the effect was not significant (z = −0.994, df =20.1, p = .332), and the predictor was dropped from the final model (likelihood-ratio test D3: p = .335).
We obtained similar results when the analysis was restricted to only complete cases (excluding samples from October). However, logistic regression with multiple imputation was generally more efficient as can be seen from the smaller confidence intervals and lower p-values (Table 6).

Also, infection load peaked in Lammefjord and Lynaes in October, which
would not have been observed in the complete cases model.  Figure 6).

| Infection load in natural host populations
Higher organic content in the sediment was associated with lower infection load, but the effect was not significant (t = −1.579, df =140, p = .117). Also, the length of the host did not have a significant effect TA B L E 5 Pooled logistic regression odds ratios for prevalence of infection

| D ISCUSS I ON
In this study, we show that ddPCR can be used as a reliable tool to quantify symbionts from whole host DNA extracts. We developed a ddPCR assay to detect and quantify an undescribed agamococcidian based on its COX1 gene, and we demonstrated its use in studying infection dynamics by documenting the prevalence of the symbiont and the infection load in four seasonally sampled populations of its host, P. elegans. ddPCR can enable further investigation of cryptic symbiotic interactions that are difficult to carry out with other molecular methods and impossible with morphological methods alone.
Here, the agamococcidian COX1 was detected in amounts ranging from 0.04 to >3000 copies per ng of total DNA extracted from hosts (mixtures of host and symbiont DNA).
As expected for ddPCR, our assay was reproducible and precise, similar to ddPCR assays used for symbiont quantification in previous studies (Koepfli et al., 2016;Srisutham et al., 2017;Wilson et al., 2015). Results were not affected by droplet generation event ( Figure 2), and variation in the concentration was small especially when the number of copies/µl was high (Figures 2 and 3). The precision of ddPCR is expected to increase with an increasing number of partitions (Quan et al., 2018). In our experiments, the droplet count varied from 10,239 to 17,553, which is typical for dropletbased methods, and it did not correlate with target concentration ( Figure 3). Additionally, the precision is affected by the average number of target molecules per droplet (λ), since an uneven distribution of template across the partitions might cause a decrease in precision for very low target concentration samples (Hindson et al., 2013;Strain et al., 2013). Our experiments showed an even distribution of targets per droplet; however, low λ and higher %CV values were detected for very low target concentrations (Tables 3 and 4).
As ddPCR results in absolute concentration and no standard curves are needed (Hindson et al., 2013), the method is particularly applicable for quantifying undescribed symbionts from total DNA extracts of hosts, when extraction of the symbionts is not feasible to provide the reference material needed when preparing standard curves (Salipante & Jerome, 2019). Importantly, our assay can detect symbionts at a very low density (0.86 copies/µl), which might realistically be expected for agamococcidians and other apicomplexans. Calculation of infection loads for agamococcidians has not been attempted previously, but Miroliubova et al. (2020)  S. pygospionis, the number of mitochondria increases with cell size (Paskerova et al., 2018). To achieve more accurate estimates of the symbiont count, a single copy nuclear marker would be better suited for the assay. Considering how well the target symbiont is known and whether genetic data from the symbiont is available, it could be difficult to design primers with sufficient specificity to use with total host DNA extracts because of conservation of gene sequences shared with the host or other related symbiont species. We used the COXI gene fragment because it was available and sufficiently different from the host COXI sequence (Supporting information 1) and because the COXI gene has shown potential for DNA barcoding of some protists (Pawlowski et al., 2012) including apicomplexans (Ogedengbe et al., 2011). It would be interesting to follow up our F I G U R E 6 Mean infection load (COX1 copies/ng total DNA) in the four populations over the sampling period. Error bars indicate the standard error of the mean Note: The references for population and month were Herslev and March, respectively. The distribution error function is normal. Because only infected hosts were used in the analysis, the sample size was 329 in total, including 5-26 individuals per population in each month. Analysis was performed without multiple imputation since the variables with missing data (length and organic content %) were not significant. study using a different assay that could provide more exact counts.

TA B L E 8 Linear regression coefficients for the log-transformed infection load
Regardless, our ddPCR assay based on COXI not only allows for detection of the undescribed agamococcidian but also provides an in-  (Kristmundsson et al., 2015). In another study, prevalence of apicomplexans infecting reef-building corals was low during summer months and increased towards fall (Kirk et al., 2013). In contrast, some marine apicomplexan infections are found to remain stable throughout seasons (Halliday-Isaac et al., 2021), suggesting that the infection dynamics of marine apicomplexans are case specific and require more thorough examination.
Apicomplexans are generally expected to be parasitic (Morrison, 2009), but how the agamococcidian studied here affects the host's fitness is not known. Aggregated distribution within hosts is typical for parasitic symbionts (Anderson & May, 1978), and we expect that the same factors influencing parasite dynamics on a general level, specifically encounter rates and host susceptibility (Schmid-Hempel, 2011), could also be relevant for the studied species. Marine apicomplexans are thought to have simple life cycles with one host species and passive oral-fecal transmission among hosts (e.g., Leander, 2007). Therefore, high host densities are expected to increase prevalence and infection loads (Anderson & May, 1978;Arneberg et al., 1998) since symbiont transmission is more efficient because of increased encounter rates (Patterson & Ruckstuhl, 2013;Rifkin et al., 2012). However, we found that when host density was highest (in May), prevalence of infection did not increase accordingly. In addition, the infection load was lower when host density was high, although the estimate was not significant. A similar non-significant negative effect of host density on prevalence was detected in apicomplexans infecting amphipods (Grunberg & Sukhdeo, 2017).
Those results, together with ours, suggest that the relationship between host density and prevalence or infection load is not necessarily straightforward and is possibly modified by other characteristics of both host and symbiont biology.
For instance, seasonality in host reproduction and reproduction strategy might affect prevalence and infection loads (Šimková et al., 2005;White et al., 1996) via changes in encounter rates and susceptibility. The studied host populations show seasonal dynamics in their reproduction, with reproductive peaks occurring in both spring and fall (Thonig et al., 2016). We observed that infection load and symbiont aggregation were high in October and November when the proportion of reproducing hosts also increases, perhaps because of an increase in food consumption to enable gamete production and consequential exposure to agamococcidian oocysts. However, we did not observe a similar peak in spring. Additionally, the host populations differ in the main mode for larval development (Thonig et al., 2016 Another host characteristic that affects symbiont encounter is host body size (e.g., Taskinen & Valtonen, 1995). Larger size is generally connected to older age or increased consumption rates, which could further increase encounter rates with the symbiont. We found that larger worms were more likely to be infected, but host size did not affect the infection load. In a system where two apicomplexans are infecting the same host species, Grunberg and Sukhdeo (2017) found that larger and older hosts had higher infection load of one apicomplexan species but not the other, leading them to propose that changes in host feeding patterns along with seasonal changes in host demographics might drive changes in symbiont abundance.
The host species in our study, P. elegans, can both filter feed and deposit feed (Anger et al., 1986), deposit feeding being more common in larger worms (Hentschel, 1998), which could explain our observation of larger worms being more likely to be infected. As the available food type can affect the grazing behavior in polychaetes (Jordana et al., 2001), filter feeding in P. elegans might be more common in spring when phytoplankton is more abundant. A change in feeding behavior with season might reflect the infection load patterns observed in this study, since hosts are more likely to encounter symbiont infective stages while deposit feeding. This could also explain why we did not observe an increase in parasite load during the spring reproductive peak. Even though reproducing hosts have higher consumption rates to support developing gametes, if the increased consumption is through filter feeding in spring, agamococcidian oocysts could be avoided. Furthermore, the symbiont's life history might also affect the encounter rates. Many marine symbionts accelerate their reproduction rate inside their hosts during warm seasons (Overstreet & Lotz, 2016). Consequently, there could be a progressive accumulation of oocysts in the sediment towards fall. Since apicomplexan oocysts are known to be very persistent (Clopton et al., 2016), oocysts produced in summer can remain viable until fall when hosts encounter rate possibly is higher because of the afore-mentioned features in host biology.
We observed differences among host populations and did not find a common seasonal pattern of infection, suggesting that environmental factors unrelated to broad seasonal changes could be influencing symbiont and host interactions. In this study, the environmental factors are all represented by the organic content % in the sediment, since other variables measured (salinity, temperature, and sediment characteristics) were correlated with host density. We chose to retain organic content in our analysis since it measures the amount of nutrients available for the host to use (e.g., Cheng et al., 1993). Higher organic content was found to lower infection load and the prevalence of infection, even though the effect was not statistically significant. Low organic content in sediment has been correlated with larger foraging area in the polychaete Streblospio benedicti (Kihslinger & Woodin, 2000) and could therefore increase the encounter rate of the host with symbiont oocysts and lead to heavier infection. The pattern could also indicate that hosts living in a habitat with higher organic content could be in a better condition and thus have better immunity against the infection. However, as the nature of the interaction between the symbiont and its host (parasitic or mutualistic) and possible immunity is not currently known, further studies are needed. In addition, it is important to remember that one of the other environmental variables correlated with organic content % might be more relevant biologically for explaining the infection patterns, and other unmeasured environmental variables cannot be discounted.

| CON CLUS ION
We have shown that ddPCR can be used effectively to detect symbionts and quantify infections from DNA extracts of hosts despite challenges of working with cryptic and undescribed symbiont species. Our assay reliably targets an agamococcidian symbiont in a marine polychaete, P. elegans, and provides specific quantification of infection load despite the fact that DNA pools could also contain DNA from other, related apicomplexan species as well as host DNA. We surveyed four natural host populations sampled seasonally and found that infection was population specific and dynamic.
A common seasonal pattern of prevalence was not detected. Since very few studies have been able to survey agamococcidians or other marine apicomplexans, our results suggest a concrete avenue towards better understanding of the interaction of marine apicomplexans and their hosts. Effective and efficient tools, such as the ddPCR assay introduced in this study, are required for monitoring the host-symbiont dynamics and the consequences of their biological interactions.

ACK N OWLED G M ENTS
This work was funded by the University of Jyväskylä Graduate School for Doctoral Studies (Department of Biological and Environmental Science) and Emil Aaltonen Foundation grant to ALH (grant number a6a412). We would like to thank Su Houyi for working on the ddPCR assay, Gita Paskerova for providing the gregarine samples, and Joel Kostensalo for helping with the modeling of missing data. We thank Hanna Hartikainen, an anonymous reviewer, and the editors for their detailed and supportive reviews and comments.

CO N FLI C T O F I NTE R E S T
Authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data and R scripts are uploaded to Dryad repository (https:// doi.org/10.5061/dryad.jwstq jq90). Sequence data generated in this study from two isolates of agamococcidian COX1 with accession numbers OK562425 and OK562426 will be available in NCBI GenBank upon acceptance.