## Abstract

The macroecological pattern known as Taylor's power law (TPL) represents the pervasive tendency of the variance in population density to increase as a power function of the mean. Despite empirical illustrations in systems ranging from viruses to vertebrates, the biological significance of this relationship continues to be debated. Here we combined collection of a unique dataset involving 11 987 amphibian hosts and 332 684 trematode parasites with experimental measurements of core epidemiological outcomes to explicitly test the contributions of hypothesized biological processes in driving aggregation. After using feasible set theory to account for mechanisms acting indirectly on aggregation and statistical constraints inherent to the data, we detected strongly consistent influences of host and parasite species identity over 7 years of sampling. Incorporation of field-based measurements of host body size, its variance and spatial heterogeneity in host density accounted for host identity effects, while experimental quantification of infection competence (and especially virulence from the 20 most common host–parasite combinations) revealed the role of species-by-environment interactions. By uniting constraint-based theory, controlled experiments and community-based field surveys, we illustrate the joint influences of biological and statistical processes on parasite aggregation and emphasize their importance for understanding population regulation and ecological stability across a range of systems, both infectious and free-living.

## 1. Introduction

For over five decades scientists have debated the driving forces underlying the macroecological scaling relationship known as Taylor's power law (TPL), which is often used to study the degree of aggregation within or among populations [1–7]. Originally outlined in a short paper to show that variance in population density increases as a power function of its mean [1], TPL has become one of the most widely verified scaling ‘rules’ in ecology, with empirical support from over 500 species of viruses, protists, arthropods, plants, birds, mammals and even humans [8–12]. Part of the appeal of TPL is its simplicity: plotted in log-log space (log*s*^{2} = log*a* + *b*log*m*), the variance (*s*^{2}) in density is a strongly linear function of the mean (*m*), for which the slope (*b*) measures the degree of aggregation [13]. Thus, higher values of *b* indicate that the variance in density (i.e. aggregation) is increasing more steeply with mean density relative to low values. Unlike other aggregation measures that change with the mean (e.g. the negative binomial parameter *k* and the variance-to-mean ratio), *b* captures variation across widely varying mean values and is thus especially suited for comparisons over time and space. While early debates focused on whether such patterns emerged from species' intrinsic behaviours or demographic stochasticity [2,3,10], subsequent extensions of TPL to a wide range of non-biological systems, including stock market fluctuations, precipitation patterns, and tornado outbreaks [12,14], have illustrated that no single explanation drives this relationship. The influence of aggregation on system stability, whether in ecological food webs or human social structures, continues to motivate research on potential applications of TPL ranging from sampling agricultural pests to protecting endangered species [15,16].

One of the most recurrent questions surrounding TPL involves the degree to which variation in the slope (*b*) reflects underlying biological processes or is instead the product of statistical or sampling artefacts [17–19]. Variation in sample size, sampling method, spatial and temporal scale, and the range of population densities considered all have the potential to influence the log-mean–log-variance relationship [2,9,20,21]. Thus, rather than comparing observations to a static null value (e.g. the slope of unity from a Poisson distribution), recent evidence highlights the importance of approaches that incorporate constraints imposed by the specific dataset [5,22]. For example, Xiao *et al*. [22] used ‘feasible set’ theory to show that the general form and slope of TPL were an emergent function of the number of individuals and groups sampled (figure 1*a*), indicating that the classical use of a Poisson null distribution is probably inadequate. Similarly, we showed in a previous study that feasible set theory could also predict the general form of aggregated host–parasite distributions [23]. Nonetheless, these previous studies show that constraint-based approaches on their own cannot typically predict either specific slope values within empirical data or capture all of the variability in an observed distribution [23–25]. The open challenge therefore becomes identifying biological mechanisms that influence aggregation (i.e. the slope of TPL) after accounting for statistical constraints imposed by the data structure.

In disease ecology, aggregation of parasites within hosts is among the most foundational and universally observed patterns [11,26–30], making infection an especially appropriate process in which to investigate biological applications of TPL [31,32]. For instance, TPL has been used effectively to characterize heterogeneity in outbreak size, superspreader distributions and the competence of reservoir hosts [30,33–35]. Building from epidemiological theory, aggregation in infection emerges from factors that cause heterogeneity in host exposure or susceptibility [29,36], including variation in host traits (e.g. immunity, sex, behaviour, genetics and age/size), spatial or temporal variation in exposure, parasite-induced mortality, and interspecific competition [28,37]. Such effects can be intensity-independent, altering the intercept of the log-mean–log-variance relationship, or intensity-dependent, altering the slope, which means that the influence varies with parasite abundance (figure 1*b*). While many studies have demonstrated the utility of TPL for characterizing parasite distributions or evaluating the influence of selected factors on aggregation [11,32,38], efforts to systematically test hypothesized drivers of such variation are often limited by data inconsistencies: meta-analytic compilations of information across divergent systems, geographical areas, sampling methods and time scales frequently introduce discrepancies. Moreover, key mechanistic processes—such as host susceptibility and parasite-induced mortality—are difficult to infer from field data alone [23,39], emphasizing the necessity of complementary experiments [11,28].

Here we integrate quantitative theory, large-scale field surveys and controlled experiments to explicitly test the roles of biological processes in driving variation in parasite aggregation. Using interactions between amphibian hosts and trematode parasites as a tractable multi-host, multi-parasite system, we measured mean infection and its variance among 11 987 individual hosts and 971 populations. After showing that feasible set theory captured both the central tendency of, and nonlinear deviations from, host–parasite TPLs, we evaluated the slope and intercept effects of empirically quantified covariates related to variation among populations, species and their environments [16,20] (figure 1). We coupled field surveys with controlled infection experiments for 20 of the most common host–parasite interactions to estimate essential epidemiological outcomes, such as host competence and parasite virulence. By comparing observed patterns in aggregation to the null model provided through feasible set theory, we show that (i) aggregation varied sharply yet consistently among host and parasite species, (ii) factors that enhanced heterogeneity within host populations (e.g. variation in body size and spatial distribution in host density) increased aggregation in an intensity-independent manner by elevating the intercept of TPL, and (iii) increases in infection competence (and especially virulence) reduced the slope of TPL, consistent with theoretical predictions but largely untested empirically. Overall, our results reveal that the shape of TPL is the joint product of both predictable, system-independent statistical properties and measurable, system-specific biological mechanisms.

## 2. Material and methods

### (a) Field sampling: community-based profiles of parasite aggregation

Between 2009 and 2015, we collected data on amphibian hosts, snail hosts and their trematode parasite infections from 173 ponds in the East Bay region of California [40,41]. We visited each wetland twice to determine community composition, the larval density of each amphibian host species, the density of snail intermediate hosts and the corresponding information on their trematode parasites (see electronic supplementary material). For each site, we measured body size (snout–vent length) and quantified infections among 10–15 late-stage larvae or recently metamorphosed individuals from all non-endangered amphibian species (*Pseudacris regilla*, *Anaxyrus boreas, Lithobates catesbeianus*, *Taricha torosa* and *T. granulosa*). We focused on infection by larval trematodes (see methods in [40]), for which four species comprise 99% of observed infections: *Ribeiroia ondatrae*, *Echinostoma* spp., *Alaria marcianae* and *Cephalogoniumus americanus*. Although each parasite uses freshwater snails as first intermediate hosts and amphibians as second intermediate hosts, they vary widely in average infection abundance, spatial occurrence among ponds, and degree of host specificity [41,42].

### (b) Experimental infections: quantifying host competence and parasite virulence

Over 4 years, we collected egg masses of each amphibian species from field sites, hatched them in the laboratory, and allowed larvae to develop to Gosner [43] stage 28 for anurans and Wong & Liversage [44] stage 2T for newts. For rough-skinned newts (*T. granulosa*), which lay solitary eggs that are difficult to find in natural habitats, we collected adults and allowed them to breed under simulated field conditions. Individual larvae were maintained in 1.5 l containers and exposed to a specific dosage of one of four trematode species: *Ribeiroia ondatrae*, *Echinostoma trivolvis*, *Alaria marcianae* and *Cephalogonimus americanus* (electronic supplementary material). The target dosages were 0 (control), 20, 40, 100, 200 and 500 cercariae per host with an average of 12 replicates per treatment, for a total of 1248 experimental hosts. At 20 days post-exposure, we recorded the fraction of surviving hosts in each treatment and measured the number of trematode meta- or mesocercariae per host. We measured competence as the percentage of administered parasites persisting after 20 days and verified that this measure was independent of dose (see also [42]). To calculate virulence, we used firth-adjusted logistic regression to estimate the effect of parasite dose on host survival over 20 days. The specific measure of virulence was the slope of the logistic regression, indicating the increase in the log-odds of mortality in response to a unit increase in initial parasite dose, when evaluated at the mean dose. Firth-adjusted logistic regression corrects for biased parameter estimates when sample size is moderate to small and the events which are being modelled are relatively rare [45].

### (c) Analysis of aggregation

In a previous study, we demonstrated the utility of applying feasible set theory to amphibian hosts and their trematode parasites, showing that the total number of parasites *P* and the total number of hosts *H* strongly constrained the shape of a host–parasite distribution [23]. Building from those results, here we used a feasible set approach to account for inherent constraints imposed by the total number of parasites and hosts sampled on the shape of TPL (figure 1*a*; see electronic supplementary material for a full description). To generate a candidate configuration in a feasible set, *P* parasites are randomly partitioned into exactly *H* hosts. One random configuration corresponds to a potential distribution of *P* parasites among *H* hosts. This random partitioning is repeated many times, and the sampled configurations are an estimate of the feasible set given *P* and *H* [22,46]. The feasible set approach makes no mechanistic assumptions but predicts the most likely aggregation level given data constraints (see electronic supplementary material), effectively recognizing that many combinations of mechanisms can generate similar levels of aggregation [25,47]. For every observed host–parasite population in the field sampling described above, we extracted the total number of parasites *P* and total number of hosts *H* and randomly sampled 1000 configurations from the feasible set. We computed the log variance of the number of parasites per host in each configuration and took the median log variance as the feasible set-predicted log variance (see electronic supplementary material for justification). We then took the difference between the observed log-variance in infection and the predicted log-variance for all observed host–parasite populations [22] (electronic supplementary material). This residual log-variance was the response variable used for all of the analyses described below. For all analyses, log refers to log_{10}.

To investigate how host, parasite and host-by-parasite interactions affected patterns of parasite aggregation, we used linear mixed effects models with an information-theoretic approach. Residual log-variance in infection (see above) was modelled as a normally distributed response and we included log-mean infection as a fixed effect with sampling year and pond identity as random intercept terms (see electronic supplementary material for justification). For a given covariate, we distinguished between two types of potential effects on the log-mean residual log-variance relationship: intercept effects and slope effects (figure 1*b*). All continuous predictor variables were scaled to have a mean of zero and a standard deviation of one. To remove the ‘forbidden zone’ [48], we only included populations that had ≥3 hosts and ≥4 parasites, which tended to eliminate host-by-parasite combinations with very low competence (260 populations removed). After further removing any populations with missing values for candidate predictor variables, our dataset consisted of 971 host-by-parasite distributions, 11 987 hosts, and 332 684 parasites. Here, each of these 971 points represents a unique combination of host species, parasite species, wetland identity and sample year.

To determine whether host and parasite species had consistent effects on the shape of TPL, we tested whether including intercept and slope fixed effects for host or parasite identity improved model fit relative to log-mean infection alone using stepwise model selection and delta Akaike's information criterion (AIC) [49]. Because this approach does not provide mechanistic insights into why particular host or parasite species affect aggregation, we subsequently substituted host and parasite identity terms for measurements of specific attributes related to host life history, characteristics of the host population, and epidemiological properties of the host-parasite interaction. At the host population level, we included log-mean host body size, heterogeneity in host body size (i.e. the residual variance in host body size after accounting for mean host body size), log-host density, and spatial heterogeneity in host density within the pond (i.e. the residual variance in amphibians captured per dipnet sweep after accounting for mean density; see [50] and electronic supplementary material for calculations). At the parasite level, we included parasite body mass, average number of cercariae released per snail, and average density of snails infected with a given parasite across sites that supported the infection (see electronic supplementary material, table S1 and [40,51,52]). We subsequently replaced host and parasite predictors with the log-transformed experimental estimates of host competence and parasite virulence for each host-by-parasite combination. We expected increasing competence or virulence to decrease the slope of TPL (figure 1) [36]. We did not include both variables simultaneously because they were positively correlated (*r* = 0.59, *p* < 0.001) and no host–parasite combination had low competence and high virulence. In the final step, we used forward model selection to re-incorporate any influential host or parasite variables and identify the best overall models (table 1; see electronic supplementary material, tables S2–S5 for full model selection). All mixed effects model analysis was performed in R (v. 3.2.4) using the ‘lme4’ package [55], and feasible set analysis was performed in Python (v. 2.7.12) using the ‘pypartitions’ package [46]. Scripts are available in the electronic supplementary material.

## 3. Results

By using a feasible set approach to account for constraints imposed by the number of sampled hosts and parasites, we were able to describe 90% of the variation in the log-mean–log-variance relationship among 332 684 parasites from 11 987 hosts (figure 2*a*,*b*). However, while the feasible set captured the central tendency of the host–parasite TPL, *P* and *H* alone were not sufficient to describe the observed patterns of parasite aggregation (21% of the observed data points fell outside of the respective 95% CIs for the feasible set-predicted log variance; figure 2*a*). Importantly, we detected strong and consistent effects of host and parasite species identity on empirical patterns of aggregation even after controlling for feasible set predictions. The TPL intercept—but not its slope—varied among amphibian species (likelihood ratio test [LRT] intercept: , *p* < 0.001; LRT slope: , *p* = 0.082; figures 2*c* and 3*a*). For instance, bullfrogs (*L*. *catesbeianus*) and newts (*Taricha granulosa* and *T. torosa*) had, on average, higher variance for a given infection intensity relative to other species (figure 3*a*). The effects of parasite identity, in contrast, manifested on both the slope and the intercept (figures 2*d* and 3*b*; LRT slope: , *p* < 0.01; LRT intercept: , *p* < 0.01), and were more important than host effects (**Δ**AIC without host: 12.6; **Δ**AIC without parasite: 139.3; compared with model 1 in table 1). Infection by the highly virulent trematode *R*. *ondatrae* was associated with the greatest reduction in slope (figure 3*b*). Including both host and parasite identity collectively accounted for 29.5% of the residual variance in aggregation not described by the feasible set and reduced AIC by 136.25 units compared with a model with only the log-mean of infection. These patterns were broadly consistent across sampling years (figure 3).

Replacing host and parasite species identity with biological traits offered insights into the potential mechanisms underlying observed patterns. Host populations with a larger average body size, more among-host variation in body size and a higher degree of spatial heterogeneity exhibited greater parasite aggregation (table 1, model 2). These were intercept effects and collectively reduced model AIC by 4 units relative to host species identity. In contrast, parasite traits had both intercept and slope effects, consistent with the effects of parasite species identity. Parasites that, on average, released more cercariae per infected snail showed an intensity-independent increase in aggregation among amphibian hosts (table 1, model 3). Moreover, increases in parasite body mass or the average density of infected snails were associated with slope-dependent decreases in aggregation (table 1, model 3).

Including traits related to the host-by-parasite interaction accounted for additional variation in parasite aggregation. Based on experimental exposures, infection competence and virulence both varied widely among the 20 tested host-parasite combinations (figure 4). *Ribeiroia ondatrae* was the most virulent parasite, yet its effects ranged from a 52% increase in *per capita* mortality risk in *A*. *boreas* to a 0.4% increase for *L*. *catesbeianus*. Incorporating these estimates into the models revealed that both terms reduced the slope of TPL (table 1, models 4–5; figure 4*a*,*b*). Host–parasite pairs with high competence (figure 4*a*) or associated with higher host mortality (figure 4*b*) both exhibited substantially lower slopes of TPL than those with low competence and low virulence. Based on the competing-model approach, the combination of host-population characteristics (body size heterogeneity and spatial heterogeneity), parasite traits (cercariae released per infected snail, parasite body mass, and infected snail density) and attributes of the host-by-parasite interaction (competence or virulence) accounted for 28% of the residual variance in parasite aggregation.

## 4. Discussion

Despite the long history and diverse applications of Taylor's power law, the identity and importance of its contributing processes continue to be debated [3,4,10,22]. In part this stems from the concurrent challenges of accounting for statistical constraints of the data [22] while also obtaining sufficiently comparable information to test alternative mechanisms simultaneously. Many studies assume that slope values greater than one—which would be the expected value according to a Poisson distribution where the mean equals the variance—are indicative of aggregation [11,34,56,57], yet the static slope provided by a Poisson is an inadequate null for determining whether this level of aggregation is ‘unusual’ given the observed data [23]. The current combination of multi-species field surveys with experimental quantification of common host–parasite interactions revealed that patterns of parasite aggregation among 971 populations were an emergent property of the host species, the parasite species, and the host-by-parasite interactions. The use of feasible set theory, which predicts the most likely aggregation patterns given data constraints, indicated that host–parasite TPLs are tightly constrained by how many hosts and parasites comprise each population [22], effectively allowing us to describe 90% of the variation in TPL shape. Unlike conventional log-mean–log-variance approaches, feasible set theory also accounts for expected non-linearities in empirical TPLs [20,22] and obviates the need for including sample size as an additional predictor [38] (figure 1*a*).

After controlling for the feasible set predictions, our results identified strong, joint effects of host and parasite species identity on infection aggregation, which collectively accounted for approximately 30% of the residual variation. These effects were consistent across a large spatial extent and over 7 years of data collection, suggesting that aggregation patterns were a function of both the species under consideration as well as its specific environment (in this case the host; see [34]). While parasite identity affected the slope and intercept of the log-mean–log-variance relationship, host species only affected the intercept (see also [10]), emphasizing the importance of testing the biological processes responsible for varying aggregation strength across different environments (i.e. host-parasite combinations). Although many studies have focused on variables that affect either the intercept or more often the slope of TPL, consideration of both intercept and slope effects on TPL provides a more conceptually integrated framework. For instance, analyses that constrain the intercept to zero (a mean infection of 1), compare log-variance estimates at a fixed value of infection or extract residuals from the log-mean–log-variance regression all preclude joint assessments of slope and intercept changes in TPL [38,57,58]. While such approaches may be necessary if the intercept is sensitive to variation in sampling method and scale, particularly when compiling data among many sources, the systematic and consistent sampling of the current study afforded an opportunity to extend the range of inference.

By linking empirical censuses of parasite aggregation with field-based measurements of population-, species- and environment-level traits, we show that characteristics of the host population—including variance in body size and spatial heterogeneity—increased parasite aggregation in an intensity-independent manner. Hosts populations with larger average body sizes exhibited greater aggregation, consistent with both the longer duration over which parasites were likely acquired and hypothesized life-history trade-offs between developmental time (or body size) and investments in immune defences [42,59–62]. Characteristics of the parasite, such as cercariae size, number of parasites released per snail host and the density of infected snails per wetland, also affected parasite aggregation. These results stemmed from the tendency of more virulent parasites to be larger in size but fewer in number [63]. Thus, low-virulence parasites, such as *Alaria marcianae*, positively affected infection heterogeneity in amphibian hosts both because its small cercariae were unlikely to induce host mortality (which reduces aggregation) and because of the inherently greater variance associated with the large number of infectious stages released per snail (figure 2*d*).

Perhaps more importantly, the current findings highlighted the influence of species interactions in determining infection heterogeneity among populations [58,64]. Previous studies have debated whether aggregation generally and the slope of TPL specifically is an intrinsic property of the species being studied or instead the product of variation in habitat, sample size and sampling methods [3,8–10]. Our results advance this debate by illustrating that, while species identity matters, aggregation is more accurately captured by variables reflective of the species-by-environment interaction. In the case of parasites, this entailed the joint roles of host and parasite species identity and the epidemiological outcomes of this interaction—namely, infection competence and virulence. Both variables had slope effects on aggregation, albeit for different reasons. Higher virulence and parasite-induced mortality likely decreased aggregation by truncating variance around mean parasite abundance, especially at higher infection loads [36,47]. These effects associated most strongly with infection by *R*. *ondatrae*, which both in the current study and past research induced mortality in a strongly dose-dependent manner [65]. For competence, host–parasite combinations with low competence were associated with larger TPL slopes relative to those with high competence (figure 4*a*). This variation is probably the product of differences in the infection process: generating a high average infection load when competence is low requires a substantial increase in parasite exposure, which will incur its own increase in variability either spatially among hosts or temporally among exposure events. Interestingly, parasite identity remained a more influential covariate than either competence or virulence (table 1), which could stem from biological processes not captured by these epidemiological predictors, such as evolutionary history, or the potential for a species-factor to better allow for nonlinear effects than a continuous variable [66].

The aggregation of parasites within hosts is considered an ‘ecological law’ in disease ecology, such that the slope of the log-mean–log-variance relationship for parasites per host often exceeds unity [26,29,30]. For instance, Shaw & Dobson [11] reported a slope of 1.55 across 263 macroparasite–host populations (s.e. = 0.037, 95% CI = [1.48, 1.62]), similar to the slope of 1.59 in the present study of 971 populations (s.e. = 0.017, 95% CI = [1.56, 1.63]) and to the mean value of approximately 1.52 reported by Xiao *et al*. [22] for free-living species (s.e. and 95% CI not available). However, efforts to identify the drivers of such variation have yielded surprisingly few generalities. While meta-analytic analyses have identified correlates associated with aggregation, including parasite specificity, coinfection, host life history, parasite developmental stage, and transmission mode [11,28,56,62,67], testing the individual and combined roles of such factors is often hindered by data availability. In a synthetic comparison of published data on aggregation by 180 parasite taxa, Poulin [38] concluded that there may be ‘little point to seeking universal causes for the remaining variation’ in parasite aggregation beyond that explained by mean infection load. As reflected by our own analyses, some of this outcome is probably methodological: the large influence of statistical processes on the perceived pattern of aggregation indicates that adequate accounting of numerical constraints via methods such as feasible set theory is often necessary before testing for biological drivers. In addition, challenges associated with a meta-analytic approach may restrict opportunities to broadly test the influence of a wide range of consistently measured mechanistic variables. For instance, in his comparative study of parasite aggregation in fish hosts, Poulin [38] was only able to include variation in host body size, parasite taxon and developmental stage as predictors; moreover, by focusing on the residuals of the log-mean–log-variance relationship without including covariate interactions with log-mean, that study was only able to test for effects on the intercept of Taylor's power law (not concurrent effects on the slope). The dataset assembled in this study is unique in the number of directly and consistently measured covariates at the population, species and species-interaction levels, with replication across multiple host and parasite species over a 7-year period, which helps explain the detection of effects not observed in previous research. Nonetheless, an important frontier will be to test these variables in other, comparably well-studied multi-host, multi-parasite systems to assess their generality and further incorporate variables such as host immunity, behavioural variation and phylogenetic history.

The ubiquity of Taylor's power law across biological and non-biological systems has led to the question of whether the shape of the log-mean–log-variance relationship contains a signature of system-specific processes beyond statistical properties of the dataset. Our analyses show that when this pattern is put in the context of the constraints imposed by the data, we can identify consistent effects of putative mechanisms on the shape of host–parasite TPLs, including both intensity-dependent and intensity-independent processes. The high predictability of parasite aggregation from constraint-based theory and the associated effects of biological mechanisms on the slope and intercept of this relationship have important implications for predicting parasite aggregation and its effects on population regulation, ecological stability and disease management [15,68,69]. While we used interactions between hosts and macroparasites as a useful and clearly delineated model system, both the results and the overall approach employed here are broadly applicable to aggregation in other biological and non-biological systems.

## Ethics

The described data were collected with the approval of the University of Colorado's Institutional Animal Care and Use Committee (protocols 1002.02 1302.01) and in accordance with sampling protocols approved by the California Department of Fish and Wildlife (SC-3683 and SC-10560), the US Fish and Wildlife Service (TE-181714), Santa Clara County Parks, East Bay Regional Parks District, East Bay Municipal Utility District, California State Parks and other local landowners.

## Data accessibility

The data and metadata associated with this article are available and have been uploaded as part of the supplementary material. They are also available on Dryad [70].

## Authors' contributions

P.T.J.J. conceived of the initial idea for the study, coordinated field and experimental studies, and organized the data; M.Q.W. and P.T.J.J. performed data analyses; and both authors contributed to the writing and editing of the manuscript.

## Competing interests

We have no competing interests.

## Funding

This work was supported by grants from the National Science Foundation (DEB-0841758, DEB-1149308 and DGE-1144085), the National Institutes of Health (R01GM109499), and the David and Lucile Packard Foundation.

## Acknowledgements

Our sincerest thanks to the many individuals who helped collect these data over the years, including especially D. Calhoun and T. McDevitt-Galles. For access to properties and logistical support during field sampling, we thank East Bay Regional Parks District, East Bay Municipal Utility District, Santa Clara County Parks, Blue Oak Ranch Reserve, California State Parks, The Nature Conservancy, Open Space Authority, Mid-peninsula Open Space and many private landowners. For discussions and feedback helpful in shaping the manuscript, we thank C. Briggs, J. Cohen, members of the Macroecology of Infectious Disease Research Coordination Network (funded by NSF DEB 131223) and three anonymous reviewers.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3861841.

- Received June 21, 2017.
- Accepted August 10, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.