Observational evidence of herbivore‐specific associational effects between neighboring conspecifics in natural, dimorphic populations of Datura wrightii

Abstract Associational effects—in which the vulnerability of a plant to herbivores is influenced by its neighbors—have been widely implicated in mediating plant–herbivore interactions. Studies of associational effects typically focus on interspecific interactions or pest–crop dynamics. However, associational effects may also be important for species with intraspecific variation in defensive traits. In this study, we observed hundreds of Datura wrightii—which exhibits dimorphism in its trichome phenotype—from over 30 dimorphic populations across California. Our aim was to determine whether a relationship existed between the trichome phenotype of neighboring conspecifics and the likelihood of being damaged by four species of herbivorous insects. We visited plants at three timepoints to assess how these effects vary both within and between growing seasons. We hypothesized that the pattern of associational effects would provide rare morphs (i.e., focal plants that are a different morph than their neighbors) with an advantage in the form of reduced herbivory, thereby contributing to the negative frequency‐dependent selection previously documented in this system. We found the best predictor of herbivory/herbivore presence on focal plants was the phenotype of the focal plant. However, we also found some important neighborhood effects. The total number of plants near a focal individual predicted the likelihood and/or magnitude of herbivory by Tupiochoris notatus, Lema daturaphila, and Manduca sexta. We also found that velvety focal plants with primarily sticky neighbors are more susceptible to infestation by Tupiochoris notatus and Lema daturaphila. This does not align with the hypothesis that associational effects at the near‐neighbor scale contribute to a rare‐morph advantage in this system. Overall, the results of our study show that the number and trichome‐morph composition of neighboring conspecifics impact interactions between D. wrightii and insect herbivores.


| INTRODUC TI ON
Plant-plant interactions have attracted the attention of plant scientists since the first suggestion that plants can communicate with one another nearly 40 years ago (Baldwin et al., 2006;Baldwin & Schultz, 1983). Both positive and negative direct interactions between plants, occurring through mechanisms such as allelopathy or volatile-mediated communication, have been described (Kalske et al., 2019;Latif et al., 2017). Indirect interactions between neighboring plants can also occur (Barbosa et al., 2009) and may be mediated by third parties in the environment such as herbivores, parasites, or natural enemies of herbivores (Barbosa et al., 2009).
Associational effects are a category of plant-plant interactions that has received much attention from plant biologists in recent years (Barbosa et al., 2009;Underwood et al., 2014). Associational effects encompass all interactions in which neighboring plants alter the likelihood that a focal plant will be attacked by an herbivore and include effects that can be positive (associational resistance) or negative (associational susceptibility) from the perspective of a focal plant.
Although associational effects are often studied in multispecies community contexts Hay, 1986), they can also occur between conspecifics within polymorphic populations (Sato & Kudoh, 2016).
While multiple factors can contribute to the strength and direction of associational effects in plant populations and communities, herbivore population dynamics are particularly influential. When herbivore populations outgrow the supply of their preferred host plants, "spillover" onto less preferred hosts can occur, producing associational susceptibility (White & Whitham, 2000). Plant apparency has also been shown to drive associational effects, as smaller plants hidden among larger neighbors have been shown to experience associational resistance (Castagneyrol et al., 2013). Associational resistance can also occur when plants induce the production of volatile compounds, which can deter herbivores from ovipositing on the induced plant and its neighbors (Zakir et al. 2013). In addition to these biotic mechanisms, abiotic factors such as microclimate variation caused by neighbors can also drive associational effects (Barbosa et al., 2009).
Datura wrightii (Solanaceae) is a system that allows for the study of interactions between conspecifics in populations that are dimorphic with respect to trichome type. This roadside weed is common in California, where two distinct trichome morphs-glandular/sticky and nonglandular/velvety-coexist in many populations (Hare & Elle, 2001). These two morphs are easy to categorize by sight and touch in the field. Previous research using a statewide sample of populations has shown that this dimorphism is maintained by negative frequency-dependent selection linked to two herbivorous insects (Goldberg et al., 2020). These two morphs do not produce significantly different herbivore-induced volatile blends (Hare, 2007), but these blends do vary based upon the herbivore species (Hare & Sun, 2011). Here, we focus on intrapopulation plant-herbivore dynamics in this system. Experimental evidence from a study of dimorphic Arabidopsis halleri revealed that associational effects occurred between glabrous and hairy individuals (Sato & Kudoh, 2016). In this system, glabrous plants lacking trichomes were more palatable to herbivores (Sato et al., 2014) and hairy plants experienced associational resistance when surrounded by glabrous ones (Sato & Kudoh, 2016). We therefore hypothesize that associational effects may exist in dimorphic populations of D. wrightii.
In this paper, we present the results of an observational study of natural D. wrightii plants across central and southern California.
Our goal was to assess the degree to which the number and trichome phenotype of neighboring conspecifics alters the likelihood of herbivory on dimorphic D. wrightii. Based upon prior work, we hypothesized that we would observe associational susceptibility (negative interactions from the perspective of a focal plant) occurring between neighboring plants with the same trichome phenotype and associational resistance (positive interactions) between neighbors with different trichome phenotypes. A scenario such as this could underlie the fitness cost associated with becoming common that occurs under negative frequency-dependent selection as it would provide locally rare morphs (e.g., those that are different from their surrounding neighbors) a fitness advantage in the form of reduced herbivory compared to the locally common morph (e.g., those that are the same morph as their neighbors).
To test this hypothesis, we measured herbivory on California D. wrightii and counted the number and phenotype of neighboring plants to look for correlations between herbivory and neighborhood size/composition (local morph frequency). We conducted field measurements three times across two years (July/August 2017; April/ May 2018; July/August 2018) so that in addition to testing for the presence of associational effects, we could also assess the degree to which it varied between and within growing seasons.

| Study system
Datura wrightii is a Solanaceous perennial shrub native to the American Southwest and Northwest Mexico. In California, these plants are dimorphic with respect to trichome type: Some plants possess nonglandular trichomes, whereas others possess glandular trichomes and feel sticky (Hare & Elle, 2001); we will refer to these as velvety and sticky, respectively. These phenotypes coexist within populations across the coastal regions of the state and differ with respect to their associated arthropod communities (Hare & Elle, 2002). Sticky plants have compromised indirect defenses (Gassmann & Hare, 2005), with the glandular phenotype conferring resistance to flea beetles (Epitrix sp.), vulnerability to mirid suckflies (Tupiochoris notatus), and proving less attractive to predatory arthropods. The trichome dimorphism is known to be governed by a single locus with classical Mendelian inheritance and is ontogenically expressed: All seedlings have glandular trichomes with the adult phenotypes exhibited following the emergence of the 5th true leaf (van . Vegetative tissues senesce at the conclusion of the growing season and re-emerge as their adult phenotype in each following year . Our study focused on interactions between D. wrightii ( Figure 1a) and five different species of herbivorous insect: Manduca sexta, Lema daturaphila, Tupiochorus notatus, Epitrix sp., and Tricobaris compacta.
All these herbivores except Epitrix are shown in Figure 1 (panels bi). Manduca sexta (Lepidoptera: Sphingidae) adults are not herbivorous but will oviposit on their solanaceous host plants-including D.
wrightii. Their larvae are voracious eaters and are capable of defoliating entire plants as the caterpillars grow from less than 1 cm to over 10 cm over the course of 2-3 weeks (Reinecke et al., 1980). Lema daturaphila (Coleoptera: Chrysomelidae) is herbivorous during both larval and adult stages and oviposit on their host plants (Kogan & Goeden, 1970). Both M. sexta and L. daturaphila are chewing insects that completely remove vegetative tissues; however, M. sexta defoliates in a regular pattern consuming all parts of the leaf; L. daturaphila on the other hand leaves irregularly shaped damage and avoids eating veins (Hare & Elle, 2002). Epitrix flea beetles (Coleoptera; Chrysomelidae) are herbivorous as both adults and larvae; however, the eggs are laid in the soil at the base of host plants and the larvae feed on root tissue (Westdal & Romanow, 1972). As such, we only assessed damage by adults, which leave numerous small puncture wounds in leaves (Hare & Elle, 2002). Tupiochorus notatus (Hemiptera: Miridae) are small piercing-sucking herbivores (as both nymphs and adults), the damage of which manifests as a distinctive yellow discoloration of leaves (Hare & Elle, 2002). Trichobaris compacta (Coleoptera: Curculionidae) adults leave small holes that are notably larger than those left by flea beetles. The larvae of this species are pith dwellers (Lee et al., 2016); thus, we did not determine their presence as it would require damaging plants and significantly disturbing their arthropod communities. Adults of all herbivore species are capable of dispersing via flight, whereas nymphs/larvae are flightless (J Goldberg, personal observation). The striking differences between the damage patterns of each herbivore allowed us to determine how much leaf area was damaged by each species individually for all focal plants included in this study. Detailed descriptions of these patterns of damage have been previously published (Elle & Hare, 2000).

| Data collection
1m of focal individual). The neighbor cutoff of 1m was selected as this distance should be sufficient for arthropods to recognize plant clusters as separate plants (i.e., plants less than 1m apart will not be recognized as individual plants; Joo et al., 2017). Only one person (JKG) measured herbivory via direct observation to avoid interobserver bias/variation. For populations with fewer than 36 plants, all plants present were sampled. The population-wide averages of these measurements were used in a prior study (Goldberg et al., 2020), but analyses of variation between individual plants in this dataset have not been previously presented.

| Statistical analysis
Our raw herbivory data contained many zeroes; therefore, we used a hurdle model approach to address the factors influencing the likelihood of herbivory being present and the intensity of herbivory on focal plants. Hurdle models are two-part generalized linear models in which the first part (the hurdle) analyzes the data using a logistic regression (positive values are collapsed to 1) and the second part uses a truncated (zeroes excluded) negative binomial regression to look at variation in the positive values only (Rose et al., 2006). For herbivore presence data (which was collected as binary), we used a binomial generalized linear mixed effect model approach. All statistics were conducted in R (R Core Team). Hurdle models were conducted using the glmmTMB() function in the glmmTMB package so that population could be included as a random variable. Population was included as a random factor in each model to account for spatial covariation.
Each dependent variable (herbivore/herbivory measures) tested was used in two models: (a) one in which the number of neighbors was included as an explanatory variable; and (b) one in which the frequency of sticky neighbors was included as an explanatory variable (for this analysis, singletons-plants with no neighbors-were excluded). Focal plant phenotype (velvety vs. sticky) was included in each model along with the interaction between it and the continuous explanatory variables. Sticky was always used as the baseline phenotype in our statistical models. In some hurdle models, damage occurrence was too low for variation in the positive values to be analyzed and some GLMMs were unable to converge with the interaction term; thus, this coefficient was excluded for analysis to proceed (noted in Tables 1-4).

| RE SULTS
Results of all models are summarized in Tables 1-4. Table 1 shows the results of hurdle models examining the relationship between herbivore damage and the total number of neighboring plants (selected results are visualized in Figures 3 and 4). Table 2 shows the results of hurdle models looking at the relationship between herbivory and the frequency of sticky plants in the immediate neighborhood (selected results shown in Figure 5). Table 3 shows the results from GLMMs looking at the relationship between arthropod occurrence and total number of neighboring plants, and Table 4 shows results of models examining the relationship between arthropod occurrence and the sticky neighbor frequency (selected results from both tables are shown in Figure 6).
Overall, focal plant phenotype was better at predicting the presence or absence of herbivore damage than was variation in neigh- coefficients were significant vs. only 2 β neighbors/sticky frequency and 3 β interaction ; Tables 1 and 2). This was also the case for truncated GLMs predicting the magnitude of herbivory (16 out of 32 β focal plant phenotype coefficients were significant vs. only 3 β neighbors/sticky frequency and 4 β interaction ; Tables 1 and 2). Focal phenotype was able to predict the likelihood or magnitude of herbivory in at least one timepoint for all herbivores included in our study (Tables 1 and 2).
We considered significant (p <0.05) β neighbors/sticky frequency or β interaction evidence of associational effects occurring at that timepoint. Consistent associational effects (occurring at all timepoints) were not observed for any herbivore, but they were observed in at least one timepoint for each herbivore except flea beetles (Epitrix sp.). In both summer 2017 and 2018, only the focal phenotype was able to predict the likelihood of L. daturaphila damage with velvety plants being more likely to be damaged (summer 2017: β focal plant phenotype = −1.439, p < 0.0001, Figure 3a, Table 1; summer 2018: β focal plant phenotype = −0.660, p = 0.00615, Figure 3c, β interaction = −0.0816 , p = 0.693; Figure 3f, Table 1).
The number of neighboring plants was found to influence the magnitude of herbivory by multiple insect species at different timepoints (Table 1; Figure 4). When all forms of herbivory were combined, we found that in the spring of 2018, sticky plants received more damage than velvety plants (β focal plant phenotype = −0.259, p = 0.00105**; The results of hurdle models comparing the likelihood and magnitude of herbivory with the frequency of sticky plants in the immediate neighborhood. Statistically significant coefficients (p < 0.05) are in bold, and nearly significant ones (0.05 < p < 0.10) are shown in italics   Figure 5f). The frequency of sticky neighbors was only able to predict the magnitude of herbivory by one herbivore:  L. daturaphila eggs (Figure 6f,g). In the summer of 2017, all variables were able to predict the likelihood of T. notatus being present on a given focal plant (β focal plant phenotype = −5.723, p < 0.001***; β neighbors = −1.531, p = 0.020*; β interaction = 2.978, p = 0.00198**, Figure 6b). We interpret this as meaning that sticky focal plants are F I G U R E 3 Results of the "hurdle" portion of hurdle models that test the ability of focal plant phenotype, total number of neighboring conspecifics, and the interaction of them to predict the likelihood of herbivory being present on focal plants. The top row (panels a, b, and c) represents the results of models with the likelihood of Lema daturaphila damage as the response variable, while the bottom row (panels d, e, and f) shows the likelihood of Manduca sexta damage as the response variable. Significant coefficients in each model are noted along with * denoting the p-value (* 0.01 < p < 0.05; ** 0.001 < p < 0.01, *** p < 0.001). Exact p-values and tests statistics are located in Table 1 0

| D ISCUSS I ON
Our results are, for the most part, consistent with previous studies of herbivore preferences for the D. wrightii morphs (Hare & Elle, 2002).
Flea beetles only rarely were found to attack sticky plants, T. notatus was consistently more likely to damage the sticky morph (and the sticky morph received more damage than velvety morphs), and M. sexta was equally likely to feed on both morphs. We found that L. daturaphila was sometimes more likely to feed on velvety plants, which was not previously observed but also not unprecedented given that F I G U R E 5 Results of analysis comparing the likelihood of a plant being damaged by any herbivore (top row, panels a/b/c) or Tupiochorus notatus (second row, panels d/e/f) and the frequency of sticky Datura wrightii plants within 1 m of focal plants. We also show the relationship between the frequency of sticky neighbors and the leaf area damaged by Trichobaris compacta (panel g). Significant coefficients in each model are noted along with * denoting the p-value (* 0.01 < p < 0.05; ** 0.001 < p < 0.01, *** p < 0.001). Exact pvalues and tests statistics are in Table 2 sticky plants have been shown to inhibit their grown rate (Hare & Elle, 2002). Our data further show that geographic heterogeneity in abundance and damage to plants is often strong for some herbivores (flea beetles, M. sexta, and T. notatus especially). We also found some signs of the local plant neighborhood as impacting the likelihood and magnitude of herbivory on our focal D. wrightii individuals.
Perhaps the most striking pattern of associational effects that we observed is that velvety plants surrounded by sticky plants appear to be more susceptible to multiple forms of herbivory.
Our data show a clear "spillover pattern" in which velvety plants in predominantly sticky patches are more likely to be attacked by T. notatus, an herbivore that primarily infests the sticky D. wrightii morph (Figures 5e and 6b). We also found that velvety D. wrightii surrounded by the sticky morph are more susceptible to infestation by L. daturaphila (Figure 6d,f,g) and are more heavily damaged by T. compacta weevils (Figure 5g). These findings-which clearly show that velvety plants do not receive a locally rare advantage-go against the hypothesis that associational effects would underlie NFDS on the D. wrightii trichome dimorphism. Indeed, these results are more consistent with the predictions of apparent competition, in which asymmetric effects of a shared natural enemy drive certain members of an assemblage extinct while allowing others to persist (Holt & Bonsall, 2017). In other words, the apparent preference of L. daturaphila for velvety plants is causing positive frequency-dependent herbivory in favor of the sticky morph. Indeed, should the observed susceptibility to herbivores by rare velvety plants result in a fitness reduction, one would expect this morph to be extirpated from predominantly sticky populations (Bonsall & Hassel, 1997). This extirpation has not been observed. Instead over the past two decades, the two morphs have continued to coexist (Goldberg et al., 2020). Taken together, our studies suggest two possibilities for why NFDS occurs over time among populations (i.e., the velvety morph increases when rare [Goldberg et al., 2020]), despite velvety plants receiving more damage when locally rare. The first is that the scale at which we measured herbivory in this study was too small to capture the full extent of herbivory on the two trichome morphs within populations in our study system. This could be compounded by variation in the density of D. wrightii (see below). The second is that the negative effect of herbivory might vary with age, being greater in small, young individuals and less in established individuals such as the focal individuals used in this study. In other words, the increase in herbivory we observed on rare velvety plants may not be enough to drive a significant reduction in fitness. This is supported by previous studies showing that D. wrightii can be exceptionally tolerant to herbivory, withstanding large portions of vegetation tissue being damaged without a reduction in seed production Hare & Elle, 2002).
The effect of neighboring plants on the susceptibility of the sticky trichome morph was more varied than the effect on velvety plants.
Our data show that the total number of neighboring plants was usually a better predictor of herbivory/herbivore presence (Figures 3b,f and 4b,e,j) than the frequency of sticky neighbors (Figure 6b). Sticky plants with more neighbors (of any phenotype) were less likely to F I G U R E 6 Results of analysis using both neighborhood metrics (total number of numbers or frequency of sticky neighbors) to predict the presence/absence of various arthropods on dimorphic Datura wrightii plants. The top panel (a) shows the relationship between total number of neighbors and the likelihood of a predatory arthropod being present on focal plants (data shown in Table 3). The remaining panels show the relationship between frequency of sticky neighbors and the likelihood of Tupiochoris notatus individuals (panels b/c), Lema daturaphila individuals (panels d/e), or Lema daturaphila eggs (panels f/g) being present on focal plants. Test statistics and p-values for panels b-g are shown in Table 4. Significant coefficients in each model are noted along with * denoting the p-value (* 0.01 < p < 0.05; ** 0.001 < p < 0.01, *** p < 0.001) be damaged by L. daturaphila (Figure 3b) or M. sexta (Figure 3f). We interpret this as an effect of herbivores spreading out across large clusters of plants and avoiding lower quality hosts (which in our system is the sticky phenotype, presumably due to noxious compounds in the exudate) in favor of higher quality host plants (velvety, in the D. wrightii system). This finding is further reinforced by the observation that M. sexta damages less leaf area on sticky plants with more neighbors (Figure 4j). This effect only appears to apply to the herbivores which infest both D. wrightii trichome morphs as T. notatus (which strongly favors the sticky morph) damaged more leaf area on plants in larger clusters (Figure 4e). Given that the predictive variable in these cases was the total number of neighbors of both trichome phenotypes, it is likely that this is a density-dependent process rather than a morph frequency-dependent process; however, because we did not strictly quantify the density of D. wrightii plants (neighborhood area in our study varies based upon the size of the focal plant), more carefully controlled observations/experimentation are required to confirm the effect of population density on herbivory in this system.
We also found evidence that likelihood of arthropod predators being found on a plant was dependent on the number of neighboring plants (Figure 6a). In this case, velvety plants were less likely to have predators on them when they had more neighbors.
This finding is interesting because it suggests that the efficacy of plant indirect defenses may be density-dependent to some extent.
Arthropod predators may avoid lower quality host plants in large clusters where higher quality options exist. Our own data show that predators are more likely to occur on sticky plants than velvety ones (Table 4, rows 16/17; not shown in any figure) and there is further evidence for this preference in the literature for other tritrophic systems (Vasconellos-Neto et al., 2007). However, this finding contradicts prior research on the D. wrightii system which showed that the sticky morph has less effective indirect defenses than the velvety morph (Gassmann & Hare, 2005). These data highlight the need for more detailed observations of the D. wrightiiassociated predator community and how predator behavior may provide asymmetric benefits to the two Datura wrightii trichome morphs in nature.
In summary, we found evidence that associational effects between neighboring conspecifics can occur within dimorphic populations of Datura wrightii. However, these effects did not match the predictions for negative frequency-dependent selection. While near-neighbor associational effects do not appear to underlie the maintenance of the D. wrightii trichome dimorphism, it is entirely possible for undetected effects to be at play. For example, herbivore populations often vary over the growing season of their hosts, and our visitation times (late April/early May for the spring 2018 visitation; late July/ early August for both summer measurements) may not correspond to the period in which these effects occur. It is also possible-given the degree to which herbivory varies from population to population (Goldberg et al., 2020)-that associational effects were occurring within some, but not all, of the populations we visited. In addition, the 1m scale at which we looked for associational effects may not match the scale upon which the processes underlying negative frequency-dependent selection are playing out and previous studies have noted the importance of scale when assessing the roles of associational effects .
Nevertheless, we showed that near-neighbor associational effects occur in the D. wrightii system, laying the groundwork for future studies into the maintenance of the balanced trichome dimorphism in California populations of D. wrightii.

ACK N OWLED G M ENTS
We would like to thank Martin and Denise Goldberg for their logistical help with the 2018 fieldwork and Jang-Dong Seo for statistical guidance. We would also like to thank Dr. Tom Pendergast and three anonymous reviewers for helpful comments on an earlier version of this manuscript. Funding was provided by the National Science Foundation (GRFP 2015195769 to JKG and DEB-1353970 to LFD) and Indiana University (Charles B. Heiser Fellowship to JKG).

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The raw data and code associated with this manuscript are available at https://doi.org/10.5061/dryad.63xsj 3v1v