Human-mediated dispersal and disturbance shape the metapopulation dynamics of a long-lived herb

long-lived Abstract. As anthropogenic impacts on the natural world escalate, there is increasing interest in the role of humans in dispersing seeds. But the consequences of this Human-Medi-ated Dispersal (HMD) on plant spatial dynamics are little studied. In this paper, we ask how secondary dispersal by HMD affects the dynamics of a natural plant metapopulation. In addi-tion to dispersal between patches, we suggest within-patch processes can be critical. To address this, we assess how variation in local population dynamics, caused by small-scale disturbances, affects metapopulation size. We created an empirically based model with stochastic population dynamics and dispersal among patches, which represented a real-world, cliff-top metapopulation of wild cabbage Brassica oleracea . We collected demographic data from multiple popula- tions by tagging plants over eight years. We assessed seed survival, and establishment and survival of seedlings in intact vegetation vs. small disturbances. We modeled primary dispersal by wind using field data and used experimental data on secondary HMD by hikers. We monitored occupancy patterns over a 14-yr period in the real metapopulation. Disturbance had large effects on local population growth rates, by increasing seedling establishment and survival. This meant that the modeled metapopulation grew in size only when the area disturbed in each patch was above 35%. In these growing metapopulations, although only 0.2% of seeds underwent HMD, this greatly enhanced metapopulation growth rates. Similarly, HMD allowed more colonizations in declining metapopulations under low disturbance, and this slo-wed the rate of decline. The real metapopulation showed patterns of varying patch occupancy over the survey years, which were related to habitat quality, but also positively to human activity along the cliffs, hinting at beneficial effects of humans. These findings illustrate that realistic changes to dispersal or demography, specifically by humans, can have fundamental effects on the viability of a species at the landscape scale.


INTRODUCTION
The two fundamental processes driving the dynamics of a plant metapopulation are dispersal among habitat patches and the demography of populations within patches. Dispersal allows the colonization of patches and can rescue declining populations (Bullock et al. 2002). As a result, dispersal limitation can cause population extinctions (Hooftman et al. 2016) or prevent colonization of new habitat (Verheyen and Hermy 2001). Consequently, any factor that impacts on dispersal is likely to have major consequences for metapopulations. Human activities are a particular contemporary cause of changes to dispersal, which are characterized as humanmediated dispersal (HMD; Wichmann et al. (2009), Bullock et al. (2018). HMD falls into two types. Indirect effects of humans (human-altered dispersal; Bullock et al. (2018)) include habitat fragmentation, which increases dispersal limitation (Lemke et al. 2015, Belinchon et al. 2017, and habitat creation, which can increase effective dispersal (Auffret et al. 2017). Direct dispersal by humans, either intentional or accidental (human-vectored dispersal; Bullock et al. 2018), may be important, particularly when humans cause long-distance dispersal (LDD), but this is rarely studied. D' Hondt et al. (2012) found evidence that dispersal of an herb via grass cuttings maintained gene flow among distant populations. Otherwise, the best evidence comes from animal metapopulations. Buk et al. (2018) showed that planned translocations of cheetahs Acinomyx jubatus helped create a viable metapopulation. In a very different example, Fountain et al. (2014) found that accidental dispersal of bed bugs Cimex lectularius by humans explained metapopulation dynamics in this species.
How the form of dispersal affects metapopulation dynamics can be examined using simulation models, which include dispersal as an explicit process. In such models of plant metapopulations, dispersal is characterized using field data (Paine et al. 2017) or by mechanistic models that represent the dispersal process (Bohrer et al. 2005, Winkler et al. 2009, Crawford et al. 2015, Castorani et al. 2017, while others represent dispersal using simple mathematical functions (Groeneveld et al. 2008). While no studies have modeled the effects of HMD compared to natural dispersal in plants, a few modeling studies have contrasted short distance dispersal and LDD, albeit on a wide variety of systems and plant types. Some studies have found that adding LDD greatly increases metapopulation size or persistence, by allowing higher colonization rates (Milden et al. 2006, Groeneveld et al. 2008). Other studies have found a variety of effects, which vary depending on model assumptions (Bohrer et al. 2005, Crawford et al. 2015. The second fundamental process, demography, can vary from patch to patch and over time. In contrast to the common assumption of the same demography everywhere in a plant metapopulation, demographic variation over both time and space can be an important driver of extinction and colonization processes (Winkler et al. 2009, Crawford et al. 2015, Halsey et al. 2016, Castorani et al. 2017. Disturbance, the removal of vegetation, is well-known as an important cause of demographic variation, often enhancing establishment of plants from algae to herbs to trees (Bullock et al. 1994, Shibata and Nakashizuka 1995, Capdevila et al. 2016. While disturbance is sometimes examined in plant metapopulation models, these generally simulate catastrophic disturbances, such as fire or logging, that kill all growing plants (although, not necessarily the seed bank) in a patch (Bohrer et al. 2005, Groeneveld et al. 2008, Meulebrouck et al. 2009, Charney and Record 2016. Smallscale disturbances, such as gap creation by grazing animals, might have important but quite different metapopulation consequences, concomitant with their positive effects on population growth. While many studies have made valuable contributions by inferring plant metapopulation dynamics from empirical patterns of colonization and extinction (Dornier et al. 2011, Belinchon et al. 2017, Manna et al. 2017, models using direct measures of processes can allow insights into the sensitivity of metapopulation dynamics to different drivers (Milden et al. 2006, Crawford et al. 2015, Halsey et al. 2016, Castorani et al. 2017. It is important to have realistic representations of dispersal and demography using field measures or validated models, and to include variability in these processes in metapopulation models (Meulebrouck et al. 2009, Winkler et al. 2009). This allows assessment of the likely role of these processes in real-world metapopulations, despite the necessary simplifications of any model. In this paper, we examine how changes in dispersal and demography affect the metapopulation dynamics of a long-lived perennial plant, wild cabbage Brassica oleracea. The system we study is a linear metapopulation growing on cliff edges (for other linear plant metapopulations see Menges (1990), Kindlmann et al. (2014), Charney and Record (2016) and Castorani et al. (2017)). As we will show, this system can be characterized as a metapopulation, having colonizations and extinctions, and an important role for dispersal among patches (Bullock et al. 2002, Freckleton andWatkinson 2002). Measuring metapopulation size in terms of both patch occupancy (Milden et al. 2006, Belinchon et al. 2017) and numbers of plants (Meulebrouck et al. 2009), we test two hypotheses: (H 1 ) human-mediated dispersal increases metapopulation size; (H 2 ) higher levels of disturbance also increase metapopulation size.

Study species and system
Wild cabbage Brassica oleracea var. oleracea (subsequently B. oleracea) is native to Western Europe where it occurs in coastal habitats, generally on calcareous cliffs (Mitchell andRichards 1979, Wichmann et al. 2008). In Great Britain, it is most abundant in the southwestern counties of Cornwall, Devon, and Dorset. It is described as a long-lived perennial, surviving up to 20 yr (Mitchell and Richards 1979). In southwest Britain, flowering is in early summer, and seeds are dispersed in late summer and begin to germinate in autumn. The small (1.5-2.6 mm diameter; Wichmann et al. 2009), smooth, spherical seeds are borne in non-explosive pods and the natural dispersal process appears to be by wind, despite the lack of obvious morphological adaptations (Mitchell andRichards 1979, Wichmann et al. 2009).
In Dorset, B. oleracea occupies a discrete region on the coast about 44 km in length, forming a narrow strip (<20 m) along the cliff edges (Fig. 1). This region is separated from other occupied regions on the south coast by tens of kilometers. Feral populations of cultivated cabbage varieties (which are not var oleracea) are found inland, but not for many tens of kilometers. Wichmann et al. (2008) found that occupancy in this region was related to local vegetation, soil properties, and aspect. They also found a strong historical signal, suggesting limited change in occupancy over the previous 70 yr.
However, B. oleracea did not appear to be occupying all potential habitat along the cliffs.

Occupancy patterns
To assess patterns and turnover in occupancy over time, we did repeated surveys through the region. We divided this coastline into 219 sections, each of 200 m length. Surveys were done during flowering in May 2003, 2006, 2007(reported in Wichmann et al. (2008 by doing a rough count of adult plants in each section (see Appendix S1 for more detail).
We examined whether human activity along the cliffs contributed to occupancy patterns, by using the InVEST recreation model. This software accesses geo-tagged photographs uploaded to the photo-sharing website Flickr, using these to calculate photo-user-days (PUD) for a location as a validated measure of the number of visitors (Wood et al. 2013). We constructed a generalized linear mixed model to assess drivers of variation in the number of years each section was occupied (0-8). We included the measures of topography, vegetation and soil nutrients used by Wichmann et al. (2008) to represent habitat quality for B. oleracea, and added the PUD for each section (values ranged from 0 to 41).

Demography
We measured seed survival and seedling establishment by sowing seed in September 2005 and again in 2006. This was done at three sites occupied by B. oleracea as well as, in 2005, three additional sites (Fig. 1). These additional sites were unoccupied by B. oleracea at the time, although they were biophysically similar to occupied sites. We originally studied these unoccupied sites to determine if B. oleracea showed changed seed and seedling demography compared to occupied sites (e.g., Belinchon et al. (2017), Dorken et al. (2017)), but this was not the case (see Appendix S2). At each site, we laid out four replicate plots and in each plot we positioned twelve 3 m transects at least 1 m apart; four disturbed and eight undisturbed (to account for likely lower establishment rates). We set up three new plots at each of the three original sites in 2006. Disturbance comprised removal of vegetation and soil to 3 cm depth and 10 cm width. This disturbance was designed to simulate the small gaps created by, for example, small mammals or death of individual plants (Bullock 2000). We sowed seeds (collected in the region the previous year) at points 10 cm apart along one-half of the transects (two disturbed and four undisturbed). We left the remaining transects unsown and used these to assess background emergence. We monitored seedling emergence and survival frequently over the next sowing. We took the mean background emergence across the four (undisturbed) or two (disturbed) unsown transects in a plot and calculated a corrected seedling emergence value by subtracting the relevant background emergence value from each sown transect in that plot.
At all six sites, we assessed seed bank survival by burying eight nylon mesh bags of 200 seeds in October 2005. We collected four bags from each site after one year and the remainder after two years. We emptied each seed bag over an individual tray of seed compost in a glasshouse and counted seedlings over three months. After this time, we applied gibberellic acid to stimulate germination and counted seedlings for a further three months. While we set up the seed bank study, we also measured the viability of fresh seeds using the same germination procedure with 24 trays of 200 seeds each.

Seed dispersal
We used the mechanistic wind and statistical humanmediated dispersal models described by Wichmann et al. (2009) for B. oleracea seed dispersal in this region. They modeled primary dispersal by wind using the stochastic simulation model developed by Soons et al. (2004) in which individual dispersal distances are generated using seed, plant, vegetation and wind characteristics. This allowed the effects of the cliffs on upward vertical wind velocities to be simulated; we used the wind dispersal model they developed that incorporated hourly extremes in vertical wind speeds. Wind dispersal deposited most (~99.9%) of seeds within 1 m of the parent, but there was a long tail with 0.01% of seeds dispersing between 50-100 m (Appendix S3: Fig. S2). Wichmann et al. (2009) also did an experimental study of human-mediated secondary dispersal by assessing the distances travelled by seeds stuck to mud on hikers' boots. Human-mediated dispersal distances were fitted well by a function that represented the rate at which seeds drop off boots as an exponential power function of distance d. This showed~80% of seeds picked up on boots fell within 15 m, with~2% travelling over 3.6 km (Appendix S3: Fig. S2). We did further work to assess the proportion of seeds that are likely to become attached to hikers' boots, the pick-up rate. We collected soil from the cliffs, placed it into 36 9 22 cm trays and scattered 250 seeds (paint-marked, see Wichmann et al. (2009)) over the soil. Wearing clean (no mud or seeds) hiking boots, a team member took a single step into the tray and we counted the number of seeds that attached to the boot. We did this 10 times for a tray and used six replicate trays, giving 60 pick-up rate estimates, which had a mean of 2.3% (range 0-16.8%).

Demographic model
We used the field data to parameterize a demographic model (Fig. 2) separately for each of the three occupied sites, Kimmeridge, Winspit, and Old Harry, in each of the six years 2001-2007, with a May-May time step. These used the adult rates specific to the relevant year and site. Because there were small differences among years and sites but strong disturbance effects (Appendix S2), we summed seedling emergence and survival values over the 2005 and 2006 studies at all six sites, but used separate values for disturbed and undisturbed treatments. We used an integral population modeling approach (Merow et al. 2014) with three stages, seeds, seedlings, and adult plants, and the rates for adult plants as function of plant size. We used linear regression of seed survival in the soil over two years to calculate initial seed viability (ϑ viable ) and annual survival (ϑ sb ), averaged among sites (Appendix S2). We calculated seedling emergence (G i ) and survival (S i ) rates from the seed sowing study for five time intervals (i = 1-5) for the 2005 study and three intervals (i = 1-3) for the 2006 study, which corresponded to sequential autumn and spring time steps over 2 yr or 1 yr, respectively (Appendix S2). Different G i and S i were used for disturbed vs. undisturbed conditions. For adult stages, following Merow et al. (2014), we calculated functions f(z) separately for each of the six annual censuses that related variables to size (z) measured as basal stem diameter (BSD) in either the current (t) or previous (t À 1) year, as relevant (details in Appendix S2). Annual survival probability was related to BSD t-1 using a second-order polynomial logarithmic function that showed a humped relationship. Proportional increase in plant size was negatively related to the log 10 (BSD t-1 ). Probability of flowering was positively related to log 10 (BSD t ). Seeds per flowering plant showed little relationship with BSD t , so we used the mean value per year.

The metapopulation model
We created a spatially explicit model in Matlab (v.7.14.0.739, MathWorks, Massachusetts) that comprised a linear metapopulation of 1,000 consecutive patches of 12 9 12 m (see Appendix S5). The 12-km extent represented a typical stretch of the study region, exhibiting varying occupancy over time (Fig. 1). As this metapopulation is somewhat continuous, we used the neighborhood size determined by seed dispersal to calculate an appropriate patch size, following Jordano (2017). Neighborhood size is 4prd % 12 m diameter, where r is the standard deviation of the natural (wind) dispersal distances (0.89) and d is the density of adults (we estimated 10 adults/m 2 ).
The model is described in detail in Appendix S3 and in Fig. 3. During each annual time step t, the population (if present) in each patch cycles through the processes described in Fig. 2, with the plants defined as seeds, seedlings of different ages and as adults in 1-mm size classes (using the mid-point rule [Merow et al. 2014]). Demographic processes are independent among patches until seeds are dispersed. Each patch of the metapopulation is assigned, at random, the demography of one of the three study sites and has this identity throughout the model run. In addition, at each year time step, one year 2001-2007 is drawn at random and all patches have their studysite-specific demography for that year. These attributes allow spatial and temporal stochasticity.
Dispersal comprises primary dispersal of all seeds by wind followed by secondary dispersal by humans of a proportion of these seeds. We used the dispersal kernel (probability density function) resulting from the wind dispersal model in Wichmann et al. (2009) to represent among patch dispersal following Colbach et al. (2001) is the distance between Cartesian x, y locations within the source patch sr and the target patch p, based on Dd = 1 m, A = Dd/(2p 9 d), the relevant proportion of the annulus for each Dd, and f(cd) is the dispersal kernel (Appendix S3). As the metapopulation is linear, this formulation reflects isotropic primary dispersal (which we assume), leading to most seeds dispersing out of the metapopulation, with only the fraction A remaining within its bounds (see Appendix S3: Fig. S3). By contrast, almost all (apart from seeds that travel beyond the ends of the metapopulation; Appendix S3: Fig. S3) secondary human-mediated dispersal stays within the metapopulation as, reflecting reality, we assume the hiking path follows the cliff top. To model secondary dispersal, we calculated (Appendix S3: Fig. S3) the proportion of wind-dispersed seeds that fall on the path of width Dd (i.e., narrower than 12 m) and combined this with the pick-up rate of 2.3%. The dispersal of these seeds to patches in either direction is where d p = (x p À x sr ) is the distance between the pickup location in the source patch x sr and the location x p in the target patch. Size is the patch size of 12 m, and Dd = 1 m. The expressionf ðdÞ ¼ b o e d b 1 , is the proportion of seeds remaining attached, following Wichmann et al. (2009), with b 0 = 0.947 and b 1 = 0.165. In each annual time step, if a patch had only a fraction of one plant or received a fraction of one seed by dispersal, that number was rounded to 0 or 1 with a probability of that fraction. This introduced a degree of demographic and dispersal stochasticity. We set a FIG. 2. The life cycle of B. oleracea with rates derived from field data (see Demographic model section). ϑ viable and ϑ sb are derived from the seed survival study, G i and S i from the seed sowing study, and the f(z) parameters from the adult tagging study. Subscript i refers to the five time intervals (i = 1-5) for the 2005 study and three intervals (i = 1-3) for the 2006 study, which corresponded to sequential autumn and spring time steps over two years or one year, respectively (Appendix S2).
Xxxxx 2020 METAPOPULATION DYNAMICS OF A HERB Article e03087; page 5 maximum of 9,000 adult plants per patch by proportionately reducing the numbers in each size class at the end of each time step, when this limit was reached. This somewhat arbitrary upper limit is roughly double the maximum local density of our field observations in 2002 (29 plants/m 2 ), noting that this observed maximum did not result in complete cover by B. oleracea in the field.

Simulations
According to the questions posed in the Introduction, we did model runs to assess (1) how changing amounts of disturbance affected the dynamics of both local populations and the metapopulation and (2) metapopulation dynamics with natural dispersal followed by human-mediated dispersal, natural dispersal only, human-mediated dispersal only, or no inter-patch dispersal (all seeds fall within the home patch). We simulated disturbance by assuming a percentage of every patch was disturbed and this proportion of seeds and seedlings underwent transitions according to the soil disturbance treatment. The remainder of seeds and seedlings had demographic values from the undisturbed treatment, and all adult plants and dispersal were not affected by soil disturbance.
Demographic model.-We ran separate demographic models for each of the three sites: Kimmeridge, Winspit, and Old Harry. The models were stochastic; for each annual time step, we drew at random one of the six years during which adult fates were measured. We drew the initial number of plants in each run from the distribution of counts in the field survey and we gave these a size distribution based on the tagged plant study for the relevant study site (Appendix S3: Fig. S1). We iterated each run for 100 time steps (yr) as a realistic time frame for population dynamics. We ran the models with no disturbance and for different amounts of disturbance up to 100%. We ran the model for each percent disturbance 1,000 times, with no upper limit to the population size. At the end of each model run, we calculated population growth rate k p .
Metapopulation model.-Each run of the metapopulation model also had 100 annual time steps. We initiated each run by populating at random 35% of the 1,000 patches, which reflected occupancy in the study region ( Fig. 1, Appendix S3), and drew the number of plants and their size distribution as for the demographic models. Each combination of percent disturbance and dispersal type was analyzed using 500 runs of the model. We calculated metapopulation growth rates in terms of the number of adults k m , and the number of patches occupied k o . In fast-growing metapopulations, occupancy of all patches and the upper limit of 9,000 plants per patch were both reached rapidly. The period of no growth after capacity was reached would lead to decreases in the calculated values of k o and k m . We avoided these artefacts in the cases where capacity was reached before the end of the run, by calculating growth rates at 95% occupancy for k o and at 50% of full capacity for k m .
We also calculated the sensitivity of the demographic and metapopulation growth rates to the individual demographic rates and the dispersal types (metapopulation only), following Hooftman et al. (2008) by calculating the proportional change in k caused by altering individual parameters by +20% and normalizing the resulting values. We focus on sensitivities for soil disturbance set at 40%, i.e., for k p, k m % 1, but also report sensitivities at other levels of disturbance in Appendix S4.

Temporal dynamics of B. oleracea on the Dorset cliffs
A majority of the 200-m sections were consistently either occupied (70 sections) or unoccupied (73) over the eight (non-contiguous) survey years. But a sizeable number of sections showed turnover, 38 being occupied in five to seven surveys, and 38 in one to three surveys ( Fig. 1; Appendix S1: Fig. S1). Forty-eight of the 76 showing turnover were not adjacent to consistently occupied sections, while 26 colonizations were of previously unoccupied sections. Furthermore, occupancy was dynamic in those sections exhibiting turnover; over the survey years we found 86 apparent colonization events and 84 extinctions (Appendix S1: Fig. S2). The measure of human activity along the cliffs, PUD, showed a positive relationship with occupancy (with both variables logged, b AE SE = 0.173 AE 0.081, P = 0.034; see Appendix 1 for the full analysis).

Modeled dynamics of local populations
The single population models indicated some differences in population growth rates among the three sites. At all sites, models parameterized for no soil disturbance generated population growth rates that indicated population decline, k p % 0.8 (Fig. 4a). Soil disturbance increased seedling emergence and survival (Appendix 2; for example, on average disturbance increased cumulative emergence over one year by 216% and survival to one year by 683%). A greater percent disturbance increased population growth rates linearly, with neutral growth (k p = 1) at roughly 35-60% disturbance depending on the site (Fig. 4a). Growth rate was most sensitive to survival of seedlings (S x ) and small adults (z ≤ 10 mm), with other elements of the life cycle having low sensitivity (Table 1). Sensitivities for populations under disturbance rates other than 40% were similar, although as disturbance, and thus k p , increased, seedling survival became more important while adult survival became less important (Appendix S4: Table S1).

Modeled metapopulation dynamics
Metapopulation growth, in terms of the number of adults in the metapopulation (k m ) was mainly determined by local population processes, whereas change in patch occupancy (k o ) reflected the dispersal mode combined with seed pressure, i.e., the number of seeds dispersing. Seeds produced in a patch were redistributed among patches, with a fraction dispersing out of the metapopulation. As we did not allow a patch to receive a fraction of one seed (fractions became 0 or 1 seed with a probability that equalled the fraction), low proportions of seeds dispersing among patches along with low seed pressure could result in no dispersal among patches.
The k m growth rate showed a similar relationship with disturbance as the single populations (Fig. 4b), with a neutral growth (k m = 1) at~35% disturbance. However, k o was neutral (k o = 1) at~10% disturbance (Fig. 4b). This is because even when disturbance levels were < 40% and so k m < 1, seed pressure from the initial metapopulation caused some increase in occupancy (Fig. 5c, d). While k m < 1 led to subsequent declines in occupancy, slow population decline under moderate disturbance levels (e.g., 15%, when k m = 0.95) meant that this initial gain was not cancelled out during the 100-yr timeframe of the model runs. (Fig. 5c). The relationships of metapopulation k values with disturbance were flatter than those of k p , reflecting the additional rate-limiting step of dispersal in growing metapopulations. Sensitivities of k m and k o to demographic processes were very similar to those for single population k p (Table 1). Again the importance of seedling survival increased, while that of adult survival decreased, with increasing disturbance (Appendix S4: Tables S2, S3). This illustrates a role for local dynamics in determining metapopulation dynamics.
Secondary dispersal by humans had a much bigger effect on metapopulation occupancy dynamics compared to primary dispersal by wind (Table 1, Fig. 5). With a 1-m path width, 0.191% of seeds underwent secondary dispersal, and 0.086% left the home patch, almost all of which fell within the metapopulation because of anisotropy. All seeds underwent primary dispersal and 0.103% left the home patch but, because wind dispersal was isotropic, only 0.010% fell into other patches. Despite accounting for a very small proportion of seeds, occupancy dynamics under only human-mediated dispersal showed almost identical patterns in occupancy increase (Fig. 5a, b) or decrease (Fig. 5c, d) as when both modes were combined. Under only primary dispersal, the increases in occupancy were much slower and the decreases more pronounced. The shorter dispersal distances under primary dispersal meant that seeds were less able to reach new patches or to rescue patches after population loss. The number of seeds dispersed by a dispersal mode had relatively low sensitivity (Table 1); showing the importance of any seeds getting to a patch far outweighed that of the amount of dispersal, since local dynamics could quickly multiply numbers in growing populations. Interestingly, metapopulation occupancy could be maintained even in the absence of dispersal among patches, because local population dynamics allowed populations to persist (Fig. 5). By contrast to occupancy patterns, in most cases dispersal mode had little effect on the numbers of adults in the metapopulation k m (Fig. 6). The pattern was slightly different under 100% disturbance, where human-mediated dispersal combined with the high population growth rates to allow rapid build-up of numbers across the metapopulation. Examination of cumulative patch colonizations and extinctions sheds more light on how dispersal mode and disturbance interacted to influence k o . Under no disturbance, colonizations caused by initial seed pressure, which were higher under human-mediated dispersal, were gradually countered by extinctions driven by negative population growth rates (Fig. 7a,b). Under 15% disturbance, the same pattern was seen, but colonizations were more and extinction rates were lower, leading to a slower decline in occupancy (Fig. 7c, d). Under higher disturbance rates, which allowed positive population growth (Fig. 7e, f), there were virtually no extinctions, and metapopulation occupancy was driven purely by the rate of colonization, as determined by dispersal mode.

DISCUSSION
Our stochastic metapopulation model of a long-lived herb was based on a wild metapopulation, and used field-derived data for demography and dispersal. The model showed strong effects of the rate of human-mediated dispersal on metapopulation size and growth, despite only a small proportion of seeds being dispersed this way, supporting H 1 . Small-scale disturbance, by enhancing seedling establishment and survival, also had a large effect on metapopulation size, supporting H 2 , ranging from decreasing metapopulations at low disturbance to rapidly increasing metapopulations at high disturbance rates. These findings illustrate that realistic changes to dispersal or demography, specifically by humans, can have fundamental effects on metapopulation processes. Indeed, we found a signal that occupancy patterns in the real metapopulation were related to human visitation rates. The model was designed to reflect the observed B. oleracea metapopulation. The field surveys revealed metapopulation signatures, with turnover in occupancy of sections between surveys, and colonization of sections that was probably by dispersal; these sections being previously unoccupied and/or not adjacent to occupied sections (and so with a low or no seedbank). This realworld situation did not resemble the classic metapopulation as first described by Levins (1969). Patches were not distinct and distributed randomly in space. Indeed, real metapopulation structures vary greatly. Many are linear, like ours; for example, following rivers (Menges 1990, Kindlmann et al. 2014, Charney and Record 2016 or coasts (Castorani et al. 2017). It can also be difficult to delineate patches in a metapopulation. One approach is to distinguish subpopulations in terms of the spatial synchrony of population processes, which requires longterm data (Cavanaugh et al. 2014). Jordano (2017) pointed out that the concept of the neighborhood size could be used to define the area of a subpopulation, based on dispersal patterns and plant density. We used this for our modeling as, given the often small dispersal distances of seeds , it ensures that not all dispersers fall within an arbitrarily defined home patch.
A related issue is determining which parts of a landscape are habitable and which are uninhabitable, the latter forming the matrix within which the metapopulation sits. Patterns of occupancy along the cliffs were related to abiotic and biotic variables, suggesting habitat quality varied. But we also found that B. oleracea seedlings established and survived following sowing in a few continuously unoccupied sections. Indeed, several studies have used seed or spore sowing, showing that unoccupied patches are indeed habitable by the species of interest (Milden et al. 2006, Belinchon et al. 2017, Dorken et al. 2017, Paine et al. 2017). This suggests it may be difficult to determine metapopulation structure simply by looking at occupancy patterns. We kept our model simple by assuming all patches were habitable, which avoided adding arbitrary amounts of uninhabitable matrix, but also reflected the occupancy patterns in shorter sections of the metapopulation (Fig. 1).
Population extinctions in metapopulations arise through demographic or environmental processes. In many metapopulation models, extinctions are forced by imposing substantial environmental changes. Catastrophic events such as fires or logging kill all growing plants and recolonization has to occur by dispersal or from the seed bank (which is "pseudo-colonization," sensu Bullock et al. 2002). In such cases, the spatiotemporal pattern of such catastrophes is a fundamental driver of metapopulation dynamics (Bohrer et al. 2005, Groeneveld et al. 2008, Charney and Record 2016, Belinchon et al. 2017, although these patterns are generally imposed by the modelers rather than derived from field data (but see Quintana-Ascencio et al. 2019). Similarly, metapopulation dynamics can be governed by succession, which decreases habitat suitability and drives extinctions, while management to earlier successional stages allows recolonizations (Meulebrouck et al. 2009).
We observed no such environmental changes, i.e., catastrophes or succession, in our study system, and extinctions in our model could happen only through demographic processes. Population growth rates fell below 1 in certain model parameterizations, leading to patch extinctions and metapopulation declines (Fig. 7). The field survey data suggested such a process, in that sections showing extinctions had many fewer plants in the previous survey (mean AE SE; 93.7 AE 12.1) than the average number of plants in other occupied sections (523.2 AE 18.6). In other parameterizations, very few populations went extinct, despite the demographic stochasticity, and the metapopulations gradually increased in size. This is not unusual; in a model of a perennial herb, Charney and Record (2016) found that demographic processes were a minor cause of population extinctions compared to catastrophic events. It is possible that demographic stochasticity for B. oleracea may be higher in the field and this would lead to more pronounced extinction and colonization dynamics. This has been achieved in other models, albeit with stochasticity levels assumed rather than drawn from field measures  (Crawford et al. 2015). Indeed, one might expect longlived perennials with seed banks to be buffered against demographic extinction, and this is apparent in metapopulation models of such plants (Milden et al. 2006, Meulebrouck et al. 2009). We investigated disturbance as an environmental driver of changes in demography. B. oleracea is described as a poor competitor (Mitchell and Richards 1979), and so it is unsurprising that local soil disturbance had a large effect on seedling establishment and survival. The individual population models and the metapopulation model that showed greatest sensitivity of growth rate (k) was to seedling survival (Table 1). This explains why the amount of disturbance had a great effect on the metapopulation growth rates, causing a switch from metapopulation decline to growth. The Dorset cliffs exhibit heterogeneity in vegetation structure, from dense cover to more open areas, which we attribute to herbivory and soil disturbance by small animals such as rabbits Oryctolagus cuniculus and voles Microtus agrestis, but also death of large B. oleracea plants. Large grazers are not present. Humans themselves also create bare areas, certainly along the path but also by off-path activities such as bird-watching or picnicking. Our disturbance experiment aimed to capture this variety by contrasting closed vegetation with the creation of small FIG. 7. Cumulative patch colonizations (red) and extinctions (blue). Along with overall occupancy (line) under different combinations of dispersal mode and percent disturbance. First column (a-e) shows secondary, human-mediated dispersal. Second column (b-f) shows primary, wind dispersal. First row (a, b) 0% disturbance. Second row (c, d) 15% disturbance. Third row (e, f) 40% disturbance. gaps. The creation of local vegetation-free disturbances, gaps, is critical for recruitment in many plant species (Bullock 2000), and through this process can increase population sizes (Bullock et al. 1994, Chandler andMcGraw 2017) and alter community structure (Myers and Harms 2009). However, we know of no other study that has included small-scale disturbance in a metapopulation model. We did so in a straightforward way, assuming every patch has the same amount of disturbance. This is simplistic, and it is possible that, in this system, there is among-patch variation in the number of small disturbances. Real disturbance dynamics can be spatially and temporally complex (Quintana-Ascencio et al. 2019), and so including any variability in our model would have required decisions, in the absence of empirical data, about the form of this variation in terms of its statistical distribution among patches, its spatial (auto)correlation, and its temporal patterns. We consider that adding a reasonable degree of spatial variation into the model would likely have increased stochasticity among runs, but would not have changed our findings. It would, however, be useful in future to explore variability in disturbance beyond the simplistic treatment in this paper and others (Bohrer et al. 2005, Purves and Dushoff 2005, Groeneveld et al. 2008, Meulebrouck et al. 2009, Charney and Record 2016. Furthermore, the very high levels of percent disturbance we modeled are unrealistic, and would likely affect survival of older plants; but these were modeled to explore this parameter space. The lower levels are more sensible, and the switch from negative to positive metapopulation growth occurs at relatively low disturbance levels;~35%. This illustrates that local, within-population processes can strongly influence metapopulation dynamics. The human vectoring of seeds had a large effect on the modeled metapopulation, specifically by increasing the rate of patch colonization. This meant that under HMD the metapopulation grew more rapidly when population growth rates were positive, and the metapopulation declined more slowly under negative population growth rates. This HMD had two aspects in contrast to natural wind dispersal; long-distance dispersal LDD and anisotropy. Based on field measures, we allowed a very low percentage of seeds to undergo HMD, but these went very long distances compared to wind dispersal. Almost all wind-dispersed seeds landed within 1 m and a tiny proportion dispersed up to 100 m. Under HMD, 2% travelled over 3.6 km. Our sensitivity analyses showed that the occurrence of these LDD events was more important than the precise numbers of seeds undergoing LDD. Metapopulation models with disjunct patches show LDD to have large effects on colonization rates and supplementing existing populations (Higgins and Cain 2002, Milden et al. 2006, Groeneveld et al. 2008. Simply by allowing a greater reach for dispersed seeds, our model shows the same applies when patches are contiguous, much as happens in models of population spread over homogeneous landscapes (Bullock et al. 2012).
HMD was particularly effective in our model because seeds disperse along the length of the metapopulation, with none lost laterally. This is because humans move along the cliff top, following the footpath. By contrast, wind disperses seeds isotropically in our model, and many seeds are lost by dispersing outside the metapopulation, to the sea on one side and farmland on the other. We assumed isotropy, but because the predominant wind directions on coasts are from sea to land and from land to sea (Greene et al. 2008), it is likely that an even greater proportion of seeds is lost laterally than in our model.
While we have shown large effects of HMD in our model, and while our earlier study (Wichmann et al. 2009) showed the process by which HMD might occur, we do not know the extent to which human activities influence dynamics in the real-world metapopulation. The relationship between human visitation rates, estimated from photo-user-days (PUD), and B. oleracea occurrence is tantalizing. We do not wish to put too much emphasis on this result as the relationship between PUD and occurrence was not particularly strong and the InVEST approach has drawbacks. There may be biases in who takes photographs and uploads them to Flickr, and where these photographs are taken (Ghermandi and Sinclair 2019). While the InVEST approach has been validated against empirical visitor data (Wood et al. 2013, Mancini et al. 2018, those analyses have used data over large areas with coarse spatial resolution, such as state parks over the United States. The reliability of PUD as an estimate of visitor numbers at the scale and resolution of our metapopulation is unclear. Allowing for these caveats, this finding suggests that humans do have a positive influence on the metapopulation, once variation in habitat quality is accounted for. The mechanism could be by HMD, but also disturbance caused by humans, or some other unknown effect of human activities. Contemporary use of the cliffs is generally for recreation and tourism, with the South West Coast Path, which traverses the metapopulation (our study area comprises roughly one-third of the Path's length in Dorset), currently having~2.7 million users per year in Dorset (The South West Research Company Ltd 2016). Human activities, and possibly their impacts on B. oleracea, reach back over centuries. Quarrying was a traditional industry here (Thomas 1998), while the cliffs have long been popular for recreation, as evidenced by their frequent mention in the works of the late 19th-century Dorset author Thomas Hardy. Humans have a variety of effects on species' spatial dynamics (Bullock et al. 2018). While most studies of native species focus on negative effects, in the case of this metapopulation, human activities seem to be beneficial.
There is growing evidence for the impacts of HMD on species' spatial dynamics (Bullock et al. 2018, Gippet et al. 2019). Our focus here has been on humans acting as dispersal vectors, and we have shown how metapopulation models can be used to explore explicitly the Article e03087; page 12 JAMES M. BULLOCK ET AL. Ecology, Vol. xx, No. xx impacts of dispersal by humans. There is a need to understand better the degree to which humans disperse different species, and its impact on populations, metapopulations and metacommunities. One way forward may be to consider dispersal by humans in the context of the Total Dispersal Kernel, TDK (Rogers et al. 2019). The TDK describes the combined contribution of all dispersal vectors to an overall dispersal kernel, and this concept encourages researchers to consider the full range of dispersal processes in a system. Ultimately, this improved understanding will only be achieved through efforts to quantify dispersal by humans in the field.