Host patch traits have scale‐dependent effects on diversity in a stickleback parasite metacommunity

1 –––––––––––––––––––––––––––––––––––––––– © 2020 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Miguel Matias Editor-in-Chief: Miguel Araújas Accepted 28 February 2020 43: 1–13, 2020 doi: 10.1111/ecog.04994 43 1–13


Introduction
A long-standing question in ecology is, 'why are some communities more diverse than others?' To address this question, biologists have sought to identify and quantify the biotic and abiotic factors that affect community diversity, and which differ across a landscape. A clear lesson from such work is that the salient factors depend on the spatial scale being considered (Cottenie 2005, Chave 2013, Meynard et al. 2013, Leibold and Chase 2017. For instance, ecosystem productivity is often negatively associated with community diversity at small spatial scales, but the trend can be reversed at large scales Leibold 2002, McBride et al. 2014). Similarly, here we show that factors affecting parasite 991 metacommunity diversity change depending on the spatial scale at which the communities are defined.
Parasites make up a large proportion of biodiversity (Poulin and Morand 2000) and have major implications for community dynamics (Lafferty et al. 2008), conservation (Thompson et al. 2010) and human health. Yet, the processes structuring variation in parasite diversity remain poorly understood (Presley 2011, Mihaljevic 2012, Moore and Borer 2012, Richgels et al. 2013, Dallas and Presley 2014, Borer et al. 2016. To better understand the rules governing the distribution of parasite richness, parasite ecologists are increasingly drawing on metacommunity theory (Grenfell and Harwood 1997, Ebert et al. 2001, Leibold et al. 2004, Seabloom et al. 2015. For parasites, each individual host represents a transient habitat patch that often supports a multispecies parasite community (an 'infracommunity' per disease ecologists). A single host population can thus be viewed as a parasite metacommunity (a 'component community'). Or, we can view each discrete host population as a habitat patch for parasites, and a collection of host populations form the metacommunity. Parasites therefore form a spatially nested metacommunity, whose diversity and composition varies among host individuals within host populations, and among host populations (Borer et al. 2016). Most research on the ecological processes governing parasite diversity have emphasized one of these spatial scales (Poulin 1997, Vidal-Martínez and Poulin 2003, Poulin 2007. There is therefore a need for more studies that bridge both spatial scales to ask whether the same processes regulate parasite diversity at both the host individual, and host population, scale? Here, we examine scale-dependent effects on parasite richness (α-diversity), using a host metapopulation of threespine stickleback Gasterosteus aculeatus on Vancouver Island. Parasite infection prevalence in stickleback is known to covary with individual hosts' diet, morphology, sex and immune genotype (Reimchen and Nosil 2001, Stutz et al. 2014, Stutz and Bolnick 2017. Parasite diversity and community composition also differ among stickleback populations (MacColl 2009, Eizaguirre et al. 2011, Poulin et al. 2011, Stutz et al. 2014, due to among-population differences in host immune genes (Matthews et al. 2010a, Eizaguirre et al. 2011, Stutz and Bolnick 2017, diet (Matthews et al. 2010a) and abiotic conditions (Simmonds and Barber 2015). These many studies of infection in stickleback yield inconsistent conclusions, perhaps because factors regulating infection have scalespecific effects. To evaluate this possibility, we considered five broadly applicable predictions concerning the host-individual and host-population 'patch' traits that structure parasite local community diversity. We focus on parasite community richness here, a companion paper (Bolnick et al. 2019) provides complementary analyses focusing on variation in the abundance of particular parasite species.

Prediction 1. Host diet has a scale-independent effect on parasite richness
In many animal species, co-occurring individuals actively prefer different prey, resources or microhabitats (Bolnick et al. 2003). Host populations also diverge in diet and resource use, often along the same ecological axes as the individual-level diet variation (Schluter andMcPhail 1992, Matthews et al. 2010a). Because many parasites are trophically transmitted (Wilson et al. 1996, Johnson et al. 2009, Cirtwill et al. 2016, host diet should correlate with parasite richness. This correlation should be consistent across spatial scales, assuming each parasite's complex life cycle is relatively invariant across the metacommunity.

Prediction 2. Host ecomorphology has a scale-independent effect on parasite richness
Diet variation often reflects variation in trophic morphology that affects prey capture and processing (Wainwright and Richard 1995). A corollary of Prediction 1 is that trophic morphology should be correlated with parasite richness. Indeed, morphology might be a more reliable covariate of parasite richness because diet measurement is noisy and often restricted to short time scales . If the diet-morphology relationship is consistent across host populations, morphology effects might therefore also act similarly across spatial scales.

Prediction 3. Host heterozygosity has a scale-independent, negative effect on parasite richness
Parasite diversity should be inversely related to the host's ability to recognize and eliminate various parasites (given equal exposure risk). Such host immune competence can be related to genetic diversity at particular loci (e.g. MHC; Wegner et al. 2003, Kalbe et al. 2009, Oliver et al. 2009), but immunity will often depend on the collective action of many loci. Thus, another measure of immunogenetic diversity might be genome-wide genetic heterozygosity (Coltman et al. 1999, Arkush et al. 2002, Šimková et al. 2008, Whitehorn et al. 2011. We predict that genome-wide heterozygosity is negatively related to parasite richness among both scales. Alternatively, parasite richness might be controlled by host immune genotypes at specific loci that affect resistance to many taxa at once. If such loci evolved in parallel across replicate populations, we should find richness associated with certain SNPs' genotypes, at both spatial scales.

Prediction 4. Host sexual dimorphism contributes only to within-population variation in parasite richness
In many species, sexes differ in feeding ecology (Shine 1991) and immunity (Zuk 1996, Nunn et al. 2009), which should influence individual-level parasite richness (Morand and Bordes 2015). Because host populations have a roughly equal sex ratio, we do not predict that dimorphism will influence among-population variation in parasite richness.

Prediction 5. Habitat contributes to among-population differences in parasite richness
Elevation, lake size and distance from the ocean are geographic features that should be experienced by all individuals living in a given site. Therefore, they should contribute only to amongpopulation differences in parasite richness (Ebert et al. 2001, Anderson and Sukhdeo 2010, Johnson and Thieltges 2010, Richgels et al. 2013, Johnson et al. 2016).

Collection
In June 2009 we collected threespine stickleback from 46 sites on Vancouver Island, British Columbia, Canada, in the historical lands of the Kwakwaka'wakw First Nations. Our sample sites included five estuaries with anadromous fish and 33 lakes and eight streams from nine watersheds (Supplementary material Appendix 1 Table A1, Fig. A1). At each site, we placed unbaited 0.5-cm gauge wire minnow traps along ~200 m of shoreline in 0.5-3 m deep water until we took 60-100 fish per site (Supplementary material Appendix 1 Table A1). Fish were euthanized in MS-222 and preserved in 10% buffered formalin after saving a fin clip in ethanol for DNA. Specimens were later rinsed and stored in 70% isopropyl alcohol after staining with Alizarin Red.

Parasite diversity
We scanned each fish under a dissection stereomicroscope to count and identify macroparasites (helminths, crustaceans, mollusks and microsporidia) to the lowest feasible taxonomic unit, following Stutz and Bolnick (2017). We focus on two parasite diversity statistics (Supplementary material Appendix 1 Fig. A2): 1) the number of uniquely identifiable parasite taxa for each fish, α ij , i.e. parasite richness in individual host i from host population j; 2) the average per-individual parasite richness, a ij , for all fish in population j.

Ecomorphology
Before dissection, we weighed all fish within 0.01 g and used digital calipers to measure standard length, body depth and body width at the pectoral fins (mm). For a random subset of 30 individuals per population, we measured three trophic traits -gape width, gill raker number and longest gill raker length -and averaged left-and right-side armor plate number. We sexed all individuals by visually inspecting gonads. Linear measurements were log transformed and size-standardized by regression on log standard length. We calculated trait means per population for population-scale analyses. Note that none of the populations exhibited a multimodal size or weight distributions that would imply a multi-year age structure, consistent with previous work on stickleback concluding that the vast majority die at the end of their first breeding season (Wootton 1984).

Diet
For a random subset of 28 populations, we analyzed fish stomach contents in the same 30 fish used for ecomorphology, recording presence of prey to the lowest feasible taxonomic level to calculate the number of observed prey taxa per fish (prey richness). Stickleback individuals and populations vary in their relative consumption of benthic invertebrates versus pelagic zooplankton (Lavin andMcPhail 1986, Ingram et al. 2011), so we use the proportion of benthic prey taxa as another diet metric. This measure was highly correlated with the first axis of a non-metric multidimensional scaling (NMDS) analysis using fish stomach contents (per the Jaccard index in the R package vegan; Supplementary material Appendix 1 Table A2). Several studies have confirmed that stomach contents are correlated with direct observations of feeding behavior, with long-term measures of diet using stable isotopes, and with trophic morphology (Matthews et al. 2010b, Snowberg et al. 2015.

Genetic diversity
With DNA from a subsample of 12 fish per population for four marine sites, 31 lakes and six streams, we used ddRADseq (Peterson et al. 2012) to obtain genotypes for 175 350 SNPs in 336 fish (mean 107 698 SNPs/individual). The protocol and bioinformatic pipeline are described in Stuart et al. (2017). We calculated genome-wide mean heterozygosity (H ij ) for each fish i in population j, and mean heterozygosity for each population ( H j ).

Statistical analyses
All statistical models were implemented in R ver. 3.5.2 (Table 1). Our primary goal in this study is to understand why host individuals and host populations differ in parasite α diversity. We fit a Poisson general linear model (GLM) evaluating how per-fish parasite richness varies as a function of host population nested within watershed (both random effects, see model M1 in Table 1). Within the majority of lakes, and in the dataset overall, a Poisson distribution was an effective description of among-individual variation in parasite richness. For each source of variation, we calculated its effect size using η 2 . Note that a companion paper (Bolnick et al. 2019) found little effect of geographic proximity (e.g. spatial autocorrelation) on the parasite community composition. So, we do not here consider spatial autocorrelation except to note that between-lake difference in mean parasite richness is uncorrelated with geographic distance as-the-crow-flies (Mantel r = 0.05, p = 0.20).

What host traits affect parasite richness variation among individual hosts?
Within each host population, we tested for correlations between individuals' phenotypic traits and their parasite richness. Due to a shared environment and ancestry, any infection differences among fish within a population should stem from individual hosts' phenotypic characteristics or genotype (i.e. Predictions 1-4). We focused exclusively on lake fish to exclude habitat differences among populations. We omitted watershed, which had 993 no effect in M1. We used a Poisson general linear mixed model (GLMM) to test whether per-fish parasite richness (α ij ) depends on a random effect of population, standard length, sex, with population random effects of length and sex (i.e. trait × population interactions; Table 1, M2). For the subset of ~30 fish with morphological trait data we ran a separate Poisson GLMM (M3) testing how α ij depends on population, sex, length, gill raker number and length, armor number, gape width and all population random effects of these traits. For the subset of lakes with diet data (Supplementary material Appendix 1 Table A1) we used a Poisson GLMM (M4) testing for effects of population, sex, length, diet NMDS1 and NMDS2, prey richness and all interactions with population (random slopes and intercepts). Full models for each of M2, M3 and M4, were evaluated against nested reduced models with AIC and likelihood ratio tests. We carried out two analyses of the effects of host individuals' genotype on individual parasite richness. First, we used genome-wide association (GWAS) to test each SNP's effect on parasite richness (M5). Specifically, we used Poisson GLMMs to evaluate how individual parasite richness depends on population, sex, length and focal SNP genotype. We retained a population-specific random effect of sex, but not length, because only the former was supported in models M2-4. We applied M5 separately for each of 39 039 SNPs (restricting our attention to SNPs scored in at least 50 fish, with minor allele frequency exceeding 10%). We applied Holm corrections to the p-values.
Second, to test effects of whole-genome genetic diversity on parasite richness, we used Poisson GLMM (M6) to relate parasite diversity to individual hosts' genome-wide mean heterozygosity, with population as a random slope and intercept, a population random sex effect and a fixed effect of length. Due to small genotype sample sizes within populations, we omit random population effects of genotype or genome-wide heterozygosity.

Does sexual dimorphism contribute to among-host variation?
The effects of sex and the sex × population interaction on parasite richness were tested in models M2-M6. To evaluate possible phenotypic mechanisms of this dimorphism, we evaluated whether diet or ecomorphology dimorphism promotes sexually dimorphic infections. First, we calculated t-statistics for sex differences in parasite richness in each lake (t[α ij ]) and t-statistics for sex differences in ecomorphology traits and diet (NMDS1 and 2, and diet breadth [prey richness]). Then we used a linear model (M7) to test for a relationship between parasite richness dimorphism and dimorphism in size, morphology and diet. We included only those traits that were, themselves, significantly sexually dimorphic and whose dimorphism varied significantly among populations. Although population constitutes the level of replication in this analysis, model M7 seeks to explain the magnitude of between-sex differences in infection (e.g. a source of withinpopulation variation).

What lake and fish traits drive parasite richness variation among host populations?
We next tested for covariates of among-population variation in population-mean per-fish parasite richness ( a i j ). logit(α ij ) ~ λ j + sex j + length j individual host sex and size M3: We focus on mean per-fish richness within populations, rather than aggregate parasite richness of a whole sample, because they are tightly correlated (r = 0.58, p < 0.0001; Supplementary material Appendix 1 Fig. A3), and the former is not biased by sample size and so does not require rarefaction. First, we tested for variation among habitats (lake, stream, estuary) while including a watershed random effect (M8). Subsequently, we focus on lakes as the level of replication (omitting stream and marine sites). We tested the effect of geographic traits on a i j (Prediction 5), using lake elevation, ocean distance and log lake size (surface area) as fixed effects (M9). Ocean distance is as-the-fish-swims distance along the river basin between the focal lake and the ocean. For the subset of lakes with ecomorphology data, we used an additional regression to test whether a i j depends on population mean length, mean gape width, mean gill raker number, mean gill raker length and mean armor plating (M10). For the subset of lakes that also have diet data, model M11 compared a i j to mean diet (% benthic prey, prey richness and NMDS1 and 2).
To test for genetic correlates of parasite richness (Prediction 3), we started with a population-scale GWAS analysis, using regression to relate a i j to the population allele frequency of a focal SNP (with watershed as a random effect, M12). Again, we used Holm corrections for multiple comparisons. We regressed a i j on population mean heterozygosity, H j (M13).
Previous analyses of microsatellites suggest that heterozygosity in these lakes depends on lake characteristics (Caldera and Bolnick 2008), so we tested a final linear model relating mean heterozygosity to lake elevation, log area and ocean distance (M14). Based on the results of models M9-14, we switched to a path analysis to account for causal relationships between predictor variables. To build the path analysis model (shown in Fig. 3), we retained significant predictor variables from the linear models, keeping at least one variable per model M9-M14.

Among-host and among-population variation in parasite richness
Per-fish parasite richness (α ij ) varied from 0 to 10 species, with a mean of 2.28 species per host (SD = 1.69, n = 4375 fish; Fig. 1). Host population explained slightly less than half the variation (η 2 = 0.426; M1: host population effect, p < 0.0001, Supplementary material Appendix 1 Table A3), with no significant watershed effect. The remaining variation arose from among-individual differences in richness. The prevalence and intensities of each parasite species, and variation in parasite community composition, are described in a separate paper (Bolnick et al. 2019).

Parasite richness covaries with individual ecomorphology
Individual hosts' parasite richness was correlated with individual traits. AIC for M2 favored retaining fish length (standard normal; Z = 11.02, p < 0.0001), sex (Z = −1.22, p = 0.22), a random effect of population, and a sex × population random slope. Larger fish tended to have higher parasite richness. Sex was not a significant main effect in M2, but there was a strong and significant sex × population interaction (p = 0.0011), i.e. the degree and direction of sexually dimorphic parasite richness varies among populations. Models with a length × population interaction were not favored by AIC model selection (Supplementary material Appendix 1 Fig. A4).
At the individual host level, we found little evidence for ecomorphology or diet effects on parasite richness. Neither gill raker length, gill raker number, nor armor plate number correlated significantly with among-individual variation in parasite richness (M3; Supplementary material Appendix 1 Table  A3). Nor were there interactions between population and ecomorphology (LRT tests p > 0.05). No diet trait affected host-level parasite richness (M4; Supplementary material Appendix 1 Table A3). There were no diet trait × population interactions. Only host size had significant effects.
Genomic analyses of individual parasite richness yielded no significant SNP-richness associations. Of 39 039 SNPs tested with model M5, 1265 were significant at p < 0.05 and 21 at p < 0.001, but no more than false discovery rate expectations. We found no effect of individual genome-wide heterozygosity on parasite richness (M6, Poisson GLMM; Z = −1.36, p = 0.174).

Among-population differences in parasite richness
Mean per-fish parasite richness ( a i j ) spanned an order of magnitude across fish populations ( Fig. 1C; M1 p < 0.001) -0.44 taxa per fish captured at Campbell River Point, to 4.72 in Gray Lake. Richness varied by habitat (M8, habitat effect χ 2 = 19.0, df = 2, p < 0.0001), being on average 2.8-fold higher in lake stickleback than in anadromous fish ( Fig. 1B; χ 2 = 14.29, df = 1, p = 0.0002) and 1.8-fold higher than stream fish (χ 2 = 7.79, df = 1, p = 0.0052). Stream stickleback had, on average, 1.5 times as many parasite taxa as marine fish (χ 2 = 3.63, df = 1, p = 0.056). Watershed had no detectable effect, so we do not consider it further. Lake geography partly influenced among-lake variation in mean parasite richness (M9). We found no support for effects of elevation (t = −0.58, p = 0.56) or lake depth (t = 0.6, p = 0.55), a weak support for an effect of log surface area (t = −1.84, p = 0.076), and a significant positive effect of ocean distance (t = 2.56, p = 0.016; Supplementary material Appendix 1 Fig. A5). A few post hoc observations are worth noting. First, although lake depth had no linear effect In (A) we use an arrow to indicate the mean per-fish parasite richness in the entire region, (B) shows the mean, 50% density and 95% density of mean within-individual richness between habitats ( a i j ), using host population as the level of replication. (C) plots a i j by population, sorted from least to most diverse. Points represent means with one standard error bars, color coded by habitat and symbols distinguishing different watersheds. on parasite richness, visual inspection of the data suggested a quadratic relationship (depth t = 4.18, p = 0.0001, depth 2 t = −3.7, p = 0.0006, Supplementary material Appendix 1 Fig. A6). Log lake surface area likewise had a quadratic effect on mean richness (area t = 2.11, p = 0.043; area 2 t = −2.64, p = 0.013, Supplementary material Appendix 1 Fig. A7). Last, although elevation had no main effect in M9, inspection of the data suggested an interaction with lake area, which was confirmed with a post hoc linear model (elevation t = 3.82, p = 0.0006, log surface area t = 3.00 p = 0.0055; elevation × area t = −3.46, p = 0.0016, Supplementary material Appendix 1 Fig. A8).
We found support for some host-population trait means explaining among-lake differences in parasite richness (M10, fish length t = 2.8, p = 0.008; size-adjusted gape width t = 2.06, p = 0.048; gill raker number t = 1.87, p = 0.072; size-adjusted gill raker length t = −0.026, p = 0.979). AIC supported a simpler model in which mean parasite richness increased with Of the 17 dimorphic cases, in 13 populations females carried more diverse infections whereas in 4 the males had more diverse infections. This variation in dimorphism direction and magnitude is associated with the direction and magnitude of dimorphism in (B) host size and (C) host diet. In each panel, dimorphism is calculated as female mean minus the male mean. So, positive values denote populations in which females were larger, ate prey scoring high on diet NMDS2 (more ceratopogonids, gammarus and stickleback eggs but fewer chironomids), and had higher parasite richness. The trends in (B) and (C) are also observed if we use Shannon-Weaver diversity of parasites (p = 0.0075 and 0.0022, respectively). mean fish length and gape width, and decreased with mean gill raker number (all typical benthic ecomorph traits). More benthic-feeding populations (higher NMDS1) tended to have more parasites per fish but this was marginally significant (p = 0.07, M11). NDMS2 and mean prey richness had no effect (p = 0.42 and 0.41 respectively).
Genome-wide association mapping (M12) revealed no significant (after Holm correction) SNP effects on host-population parasite richness. Mean genome-wide heterozygosity had no detectable association with population mean parasite richness (M13, t = 0.925, p = 0.364).
The variables considered above are likely to be inter-related. For example, heterozygosity increases with lake area (Caldera and Bolnick 2008) but decreases with distance from the ocean and elevation (M14; p = 0.029, 0.012, <0.001 respectively). To account for this covariance among predictors, we used path analysis to partition direct and indirect effects (Fig. 3, 4, Supplementary material Appendix 1 Table A4), taking the significant effects from the a i j -based linear regressions above. The path analysis explained 45.6% of the variation in mean per-fish parasite richness among lakes. Lake area had no significant direct effect on parasite richness (r = 0.053, p = 0.822), but had indirect effects via fish heterozygosity and diet. Specifically, larger lakes have more genetically diverse fish (r = 0.366, p = 0.011), and heterozygosity has a positive effect on parasite richness (r = 0.360, p = 0.038). Stickleback in larger lakes also consumed a smaller proportion of benthic prey (r = −0.777, p < 0.001), which conferred higher parasite richness (r = 0.570, p = 0.013). The indirect negative effect of lake size via fish diet exceeded its positive effect via fish heterozygosity (r = −0.311 versus 0.130). In simple bivariate correlation tests, lake distance from ocean had a positive effect on parasite richness (r = 0.310, p = 0.017). However, this positive correlation is mediated via an indirect positive effect of ocean distance on heterozygosity (Supplementary material Appendix 1 Table A4; r = 0.553, p = 0.002), which has a positive partial correlation with parasite richness. Lakes farther from the ocean also had larger fish (r = 0.52, p = 0.001), but in this analysis, mean per-lake fish size had no further relevance to fish diet (r = 0.16, p = 0.111) or parasite richness (r = 0.075, p = 0.635). Last, higher elevation lakes had less heterozygous fish (r = −0.74, p < 0.001), indirectly reducing richness.

Discussion
Parasites form diverse multi-species assemblages in nested metacommunities across individual hosts, and among host populations (Seabloom et al. 2015). Key questions in parasite ecology, and ecology more broadly, are 1) what processes dictate this metacommunity diversity, and 2) are these processes scale-dependent? Here, we present a case study using macroparasites from a stickleback metapopulation to show that parasite richness varies widely among host individuals and across host populations. In both cases, the variation in richness spans an order of magnitude. We have identified ecological factors contributing to this variation in parasite diversity and have shown that mostly different factors act at small and large spatial scales (Fig. 5).

Individual scale: why does parasite richness differ among host individuals?
Individual stickleback hosts within a population experience approximately similar abiotic conditions and biotic prey communities. Thus, variation in parasite richness among host individuals is either stochastic or due to differences in individual traits. We confirmed that stickleback individual size, gape width, diet and sex were correlated with individual parasite richness, in sometimes subtle ways, but did not find any effect of genotype or heterozygosity.
Larger individual stickleback carried more parasite taxa, as in other fish species (Calhoun et al. 2018). Larger fish may be older, having more time to accrue infections or having senescing immune systems (Zelmer 2014). Or, larger stickleback might be capable of eating more prey due to larger stomachs, increasing overall intake rates of parasites. A third possibility is that larger fish eat particular kinds of prey that promote parasite diversity. For instance, in stickleback and many other fish, larger individuals are more benthivorous. Although we do not have an a priori expectation that benthic-feeding fish would have greater parasite richness, we do also find that fish with larger gapes (a typical benthic trait) also have higher parasite diversity. Moreover, our sexual dimorphism analysis suggests that parasite richness is higher in whichever sex has the more benthic diet and is larger. Thus, we conclude that individual size and trophic ecology affect individual hosts' parasite richness. This confirms that habitat patch characteristics (i.e. individual host traits) modify community assembly, likely by modulating the rate of colonization.
Diet effects are consistent with previous studies in other species that found that individual diet influences 998 parasite infection (Amundsen et al. 2003, Stutz et al. 2014, Cirtwill et al. 2016, Hayward et al. 2017, and with the general notion that parasite communities are structured by colonization dynamics (Worthen andRohde 1996, Zelmer 2014).
Males and females differ in diet and microhabitat use in most species (Shine 1991, McGee and Wainwright 2013, Reimchen et al. 2016. Thus, individual encounter rates with different parasites should vary by sex (Reimchen and Nosil 2001). In stickleback, we found that sex-biased parasite richness covaries with sexual dimorphism in ecomorphology traits, implying a role for encounter filters. Across all lakes, females tended to be larger than males, consume more large prey, and have higher parasite richness. In lakes with larger dimorphism in body size or diet, parasite richness was also more dimorphic. The exceptions support our conclusions: in a few lakes where diet or size dimorphism is reversed, parasite richness dimorphism is also reversed. Sexual dimorphism in immune function (Zuk 1996, Nunn et al. 2009), which we did not measure, could contribute to the residual variation in sex dimorphism of parasite richness.
We found no single-locus or genome-wide genetic control of per-fish parasite richness (α ij ). This is perhaps because parasite richness is an aggregate measure arising from many parasite species' interactions with the host, which might be controlled by separate gene(s). In the companion study (Bolnick et al. 2019), we report many genetic associations with infection by each common parasite. But, because these associations tend not to involve the same genes, there is no overall genetic effect on parasite richness.

Population scale: why does parasite richness differ among host populations?
Parasite richness varied substantially among host stickleback populations with habitat, geography, ecomorphology and genetic diversity. The strongest effect was habitat: parasite richness was nearly three times higher for lake stickleback Mean proportion benthic prey Mean parasite richness Figure 4. Scatterplots of the correlations tested in the path analysis illustrated in Fig. 3, with linear regression fit lines and confidence intervals for the raw data. Each panel is labeled (a) through (m) to match the path arrows in Fig. 3, omitting path l which was not significant. Above each panel we provide the partial correlation and its significance, from the path analysis (Supplementary material Appendix 1 Table  A4 for details). In cases where the simple bivariate correlation yields a different trend than the path analysis (e.g. when one is significant and the other is not), we provide both the bivariate correlation then the path results. than their anadromous relatives; stream fish were intermediate, consistent with previous studies (MacColl and Chapman 2010, Eizaguirre et al. 2011. Given that stickleback are relatively recent (post-glacial) colonists of freshwater, this higher infection rate indicates that the colonization process increased infection risks (Weber et al. 2017), rather than providing a means of escape from historical enemies (Grunberg et al. 2019). It remains unclear to what extent these habitat differences are a function of abiotic conditions, availability of intermediate or terminal hosts, or habitat differences in dilution effects (e.g. higher non-host fish species richness in the ocean; Becker et al. 2014).
For lake populations, lake geography affected parasite richness, which was higher in mid-sized lakes (with intermediate depths), farther from the ocean, and at higher elevation. Some geographic effects acted indirectly via host population traits. For example, stickleback in larger lakes have more limnetic diets on average (Lavin and McPhail 1986), and this diet shift is associated with reduced parasite richness; i.e. lake size had a net negative indirect effect on parasite richness via diet. At the time we sampled parasites, average water temperature in the top 10 m of water was not correlated with lake size, depth and elevation, but we cannot rule out that more generally water temperature might be a causal link between lake geography and parasite community diversity.
This negative effect was lessened somewhat by a positive but weaker indirect effect of lake size through host heterozygosity. Larger lakes support more genetically diverse stickleback populations (Caldera and Bolnick 2008). This increase in mean heterozygosity was associated with higher mean per-fish parasite richness, inconsistent with the oftcited immunological benefit of genetic diversity , Joly et al. 2008, Kaunisto et al. 2013. This surprising effect might be explained if we consider a reversed causeeffect relationship: richer parasite communities might select for greater genetic diversity. But, such selection is unlikely to affect genome-wide heterozygosity, because 1) most SNPs should be effectively neutral, and 2) heterozygosity is mostly a function of lake geography.
Heterozygosity was also influenced by lake elevation and distance to ocean. Higher-elevation lakes tend to have lower heterozygosity, consistent with expectations of a stronger bottleneck during colonization. Path analysis suggests that elevation reduces parasite richness indirectly via reduced heterozygosity, but there was no direct elevation-richness Figure 5. Summary of inferences at the among-host and among-population scales, indicating statistically significant effects (with checkmarks), non-significant effects (NS) and effects not operating at a given scale (NA). * Sex differences in diet contribute to between-sex differences in parasite diversity within populations. correlation. Likewise, lake distance from the ocean only had an indirect effect on richness via heterozygosity.
Our results do not conform to general predictions from island biogeography and basic metacommunity theory. All else equal, communities should be more diverse on habitat patches that are larger, or closer to other such patches (MacArthur and Wilson 1963), as demonstrated for human parasites on islands in Jean et al. (2016). Larger lakes are akin to larger islands, but in our study, lake size had no net relationship to parasite richness due to conflicting effects via host diet and heterozygosity. Furthermore, we find that more isolated stickleback populations with longer distances to the ocean have richer, rather than poorer, parasite communities. We hypothesize that these unexpected effects could be driven by host evolutionary genetics.

Scale-dependent and -independent factors affecting parasite richness
To what extent do the variables structuring parasite diversity generalize across scales (Fig. 5)? Only two traits showed similar effects across scales. Larger individual fish had richer parasite infections, and populations with larger average fish had higher mean parasite richness. We also detected a marginally significant positive effect of gape width on parasite richness at both scales, which suggests the trend at each scale is real rather than type I error. These effects of size and gape confirm Prediction 1.
Direct measures of diet had inconsistent effects across spatial scales. At the larger spatial scale, more benthivorous (lower NDMS1) stickleback populations had richer average per-fish parasite communities. No such trend was found at the individual scale. Diet sexual dimorphism (acting within populations) did correlate with parasite richness dimorphism, but this involved NDMS2. Thus, diet does influence parasite richness, but does so inconsistently across scales. This does not support Prediction 2.
A scale-specific relationship was also seen between genomic heterozygosity and parasite richness. Some studies suggest that heterozygosity might increase a host's repertoire of immunological tools, thereby lowering richness (Coltman et al. 1999, Arkush et al. 2002, Whitehorn et al. 2011. Indeed, fish species with higher heterozygosity had fewer parasite species , Joly et al. 2008, Kaunisto et al. 2013. Our data did not support such a trend at the individual-level, despite high statistical power, refuting Prediction 3. The lack of an individual-level effect undermines the notion that heterozygosity acts through individual immunogenetic diversity. We did, however, find a population-level effect of heterozygosity, though in an unexpected direction (parasite richness increased with host population mean heterozygosity).
So, why does heterozygosity matter at the population scale? One possibility is that larger host populations can both support more diverse parasite communities and maintain higher heterozygosity. Parasite diversity can be related to host population size, following a species-area relationship (Bagge et al. 2004, Zelmer 2014. However, our path analysis found no support for a direct effect of lake size on parasite richness. An alternative explanation is that fish species exposed to more diverse parasites might be subject to stronger balancing selection and evolve higher heterozygosity (Hamilton 1982, Berenos et al. 2010. But, this adaptationist interpretation should not impact genome-wide heterozygosity, only those loci involved in immune defense and sites linked to them.

Conclusions
Because parasites are such a large component of biological diversity and have large effects on the communities in which they are embedded (Lafferty et al. 2008), the field of ecology needs to understand the biological processes regulating the distribution, abundance and diversity of parasites. Here, we present a case study illustrating the highly multivariate and scale-dependent nature of these processes. A companion paper shows similarly scale-dependent and multivariate effects on the prevalence, and covariance, between particular parasite species that constitute the communities examined here (Bolnick et al. 2019). The macroparasite community of threespine stickleback is structured by among-individual variation in host sex, size, morphology and diet, and amongpopulation differences in host size, morphology, diet, genetic diversity and habitat. Some of these variables (size, gape width) act consistently across individual-and populationscales. In general, it seems that more benthivorous stickleback have richer parasite communities, a trend observed between sexes, and among populations. Other variables are scale-dependent and contribute to parasite differences only among individual hosts (sex) or only among host populations (heterozygosity, habitat). Such scale-dependent effects are expected in metacommunities of all kinds (Leibold et al. 2004). This scale-dependence may help explain inconsistent results among studies conducted at disparate scales. The implication of these findings is that studies of species distributions and abundances need to draw scale-specific inferences.
Permits -Collection and animal handling were approved by the University of Texas IACUC (07-032201) and a Scientific Fish Collection Permit from the Ministry of the Environment of British Columbia (NA07-32612).