Reducing herbivory in mixed planting by genomic prediction of neighbor effects in the field

Domain: earth_science
Published: 2024-10-07
Editorial outcome: minor_revision

Abstract
Genetically diverse populations can increase plant resistance to natural enemies. Yet, beneficial genotype pairs remain elusive due to the occurrence of positive or negative effects of mixed planting on plant resistance, respectively called associational resistance or susceptibility. Here, we identify key genotype pairs responsible for associational resistance to herbivory using the genome-wide polymorphism data of the plant species Arabidopsis thaliana. To quantify neighbor interactions among 199 genotypes grown in a randomized block design, we employ a genome-wide association method named “Neighbor GWAS” and genomic prediction inspired by the Ising model of magnetics. These analyses predict that 823 of the 19,701 candidate pairs can reduce herbivory in mixed planting. We planted three pairs with the predicted effects in mixtures and monocultures, and detected 18–30% reductions in herbivore damage in the mixed planting treatment. Our study shows the power of genomic prediction to assemble genotype mixtures with positive biodiversity effects.

Introduction
Genetic diversity is increasingly recognized as a critical facet of biodiversity that should be conserved as a provider of various ecosystem services. In terrestrial ecosystems, plant genotypic diversity can increase plant resistance to natural enemies as the number of plant genotypes in a contiguous group of plants, namely a stand, increases. However, such a stand of multiple plant genotypes does not always result in positive outcomes. Identifying beneficial pairs from a mixture of genotypes helps us design a desirable mixture that improves stand-level properties, such as resistance to herbivores or pathogens.

In anti-herbivore defense, both positive and negative effects of mixed planting on stand-level resistance have been reported in the literature. These phenomena are referred to as associational resistance and associationalsusceptibility, respectively. Associational resistance and susceptibility involve ecological interactions that are mutually non-exclusive, such as plant–herbivore, plant–plant, and plant–carnivore interactions. In plant-herbivore interactions, associational resistance or susceptibility occurs when chemical and physical plant traits jointly repel or attract herbivores to neighboring plants. For instance, plant odor and apparency may affect the settlement of herbivores, while physical barriers and toxic metabolites may alter herbivore behavior and growth after settlement. Plant–plant interactions may modulate the expression of these plant traits through volatile-mediated communications or direct competition. Plant–carnivore interactions may also lead to associational resistance or susceptibility when a mixture of plants compared with the component monocultures attracts more or fewer natural enemies of herbivores. The complexity of associational resistance and susceptibility makes it difficult to distinguish between positively and negatively interacting genotype pairs for anti-herbivore resistance. Nonetheless, the application of associational resistance to agriculture is becoming more and more important to reduce the use of insecticides.

Despite long-standing interest, few studies have employed genome-wide polymorphism data to investigate stand-level properties in biodiversity research. Other than anti-herbivore resistance, some pioneering studies have used genome-wide association studies (GWASs) to dissect the genetic basis underlying stand-level growth within the model plant Arabidopsis thaliana. For example, studies on 98 A. thaliana genotypes have reported quantitative trait loci that mediate neighbor genotype effects on plant growth and consequently mitigate belowground competition. However, these studies require considerable effort in pairwise cultivation to sufficiently control neighbor composition among many genotypes. Due to this combinatorial cultivation effort, previous studies focused on limited pairs between ten focal and 98 counterpart genotypes in a controlled environment. These practical limitations remain an obstacle to enabling large-scale GWAS and thereby identifying the most beneficial pairs in anti-herbivore resistance in field environments.

In this study, we aimed to predict key genotype pairs that reduce herbivory by combining genome-wide single nucleotide polymorphisms (SNPs) in A. thaliana with a recently developed GWAS method named “Neighbor GWAS” (Fig. 1). This GWAS method has the same structure as the Ising model of magnetics, which has been applied to various spatial patterns in biology such as cellular development, animal skin colors, and forest gap dynamics. Neighbor GWAS employs this widely applicable model for spatial variation in quantitative trait to distinguish locus-wise positive or negative interactions between focal and neighbor individuals. Neighbor genotype effects estimated by the Neighbor GWAS method determine whether a mixture of two genotypes at a given locus alter phenotypic values at population level, thereby distinguishing SNPs with positive or negative effects. The practical advantage of Neighbor GWAS lies in its applicability to randomized mixtures of many genotypes, providing a suitable method to analyze how neighbor genotypes shape herbivore damage across space. To apply this method, we first planted replicated individuals of 199 A. thaliana genotypes at two field sites and observed herbivore damage and naturally emerging communities of herbivores and associated insect species, which were analyzed as extended phenotypes of the plants. We then used Neighbor GWAS as a tool to quantify the phenotypic variation explained by neighbor genotype effects and to conduct GWAS of neighbor genotype effects on herbivore damage and insect community composition. Genome-wide regression incorporating neighbor genotypes was used for genomic prediction of associational resistance or susceptibility out of all possible 19,701 pairs among the 199 genotypes. To test associational resistance, we finally planted three prospective beneficial pairs in mixtures as well as in monocultures. Our joint study using the recently developed method and large-scale field experiments identified key genotype pairs responsible for associational anti-herbivore resistance.

Results
Herbivore damage on field-grown Arabidopsis thaliana
To enable GWAS of herbivore damage, we planted A. thaliana genotypes in a randomized block design in two experimental gardens over two years. This allowed us to monitor the individual number of 18 insect species on approximately 6400 individual plants (≈199 genotypes × 8 blocks × 2 sites × 2 years) at a native (Zurich, Switzerland) or exotic (Otsu, Japan) field site. We quantified herbivore damage as the number of leaf holes in Zurich and leaf area loss in Otsu because the major herbivores in Zurich were flea beetles, and those in Otsu were caterpillars of diamondback moths or small white butterflies. Typically multiple insect species were found on each individual plant (average = 2.1 and 1.8 in Zurich and Otsu, respectively). The fact that the major herbivore species observed in our study were specialists of Brassicaceae led us to classify these insects into external chewers (i.e., herbivores that harbor outside and chew plant tissues) and other herbivores. The former is known to trigger the defense pathways mediated primarily by the plant hormone jasmonic acid, and the latter includes species that induce defense through salicylic acid. To identify the insect functional groups responsible for the herbivore damage, we quantified three extended phenotypes of insect communities: the individual number of external chewers (flea beetles in Zurich or caterpillars in Otsu), the individual number of other herbivores (aphids, thrips, stink bugs, and leaf miners), and the total number of insect species per individual plant. All four phenotypes exhibited quantitative phenotypic variation among the individual plants, providing suitable target phenotypes for GWAS.

Before using the Neighbor GWAS, we performed a standard GWAS to examine focal genotype effects on herbivore damage and the three extended phenotypes of insect community composition. For all four phenotypes, we found significant heritability among plant genotypes at both sites (likelihood ratio test, χ2 > 19.4, d.f. = 1, p < 0.01). Our previous study detected single-gene effects of the trichome developmental gene GLABRA1 (GL1) on resistance to herbivore damage made by flea beetles; thus, two glabrous mutants were included to test this effect. As expected, we detected significant effects of the mutation in the GL1 gene on the herbivore damage in Zurich (p < 0.05 after Bonferroni correction of multiple testing for SNPs). Although other studies have reported significant effects of the glucosinolate genes GS-OH and MAM1 on herbivory, none of the measured phenotypes showed significant peaks near these glucosinolate genes. Presumably, this was because the majority of herbivores observed in this study were specialists and thus overcame the glucosinolate defense. The results of the standard GWAS agreed with previous evidence for physical defense, whereas the herbivore damage observed in our study was unlikely to be attributable to the known variation in defense due to glucosinolates.

Genome-wide neighbor effects contributed to shaping herbivore damage
To test whether the observed herbivore damage was affected by neighbor genotypes, we quantified the phenotypic variation explained (PVE) by neighbor genotypes using the Neighbor GWAS model that considered neighbor genotype effects besides the focal genotype effects. This PVE corresponds to heritability attributable to neighbor genotypic effects as well as focal genotype effects. Compared with focal genotype effects alone, we found that the inclusion of neighbor genotypes explained a significant fraction of the phenotypic variation in the observed herbivore damage (likelihood ratio test, χ2 > 7.4, d.f. = 1, p < 0.01). Permutation tests confirmed that the observed fraction of herbivore damage explained by neighbor genotype effects was greater than that of randomly permuted neighboring plants in Zurich (permutation test, p < 0.05). This was also supported by a permutation scheme called genome rotation that randomized population genetic structure (permutation test, p < 0.05). These results show the importance of neighbor genotypes in shaping herbivore damage in the field.

Similarly, we quantified the phenotypic variation explained by neighbor genotypes for the three extended phenotypes of the insect community composition to examine which types of observed herbivores were the most influenced by neighbor genotypes. Specifically, we asked whether neighbor genotypes were more likely to influence insect species with higher mobility. The vast majority of external chewers in Zurich were flea beetles that could jump between plants, and the individual number of external chewers on focal plants was significantly influenced by neighbor genotypes (likelihood ratio test, χ2 = 6.0, d.f. = 1, p < 0.05). In contrast, the contribution of neighbor genotypes to the individual number of external chewers on focal plants was not significant in Otsu, where the major external chewers were sedentary caterpillars that could not move quickly between plants. In Otsu, instead, flower thrips that could move between flowering stems were abundant in the category of other herbivores, and accordingly, the individual number of other herbivores was significantly influenced by neighbor genotypes (likelihood ratio test, χ2 = 5.6, d.f. = 1, p < 0.05). Consistent with the significant influence of neighbor genotypes on either external chewers in Zurich or other herbivores in Otsu, the total number of insect species at both sites was affected by neighbor genotypes (likelihood ratio test, χ2 > 4.8, d.f. = 1, p < 0.05). These patterns of insect communities suggest that neighbor genotypes are more likely to influence mobile herbivores than sedentary ones.

We then conducted GWAS of neighbor genotype effects on the herbivore damage and the three extended phenotypes of insect communities. To attribute phenotypic variation to each SNP, we mapped the statistical significance of the neighbor genotype effect on the A. thaliana genome and depicted Manhattan plots. This GWAS did not detect any significant SNPs for any of the four phenotypes at each site (p > 0.1 after Bonferroni correction of multiple testing for SNPs). This was also supported by permutation tests that shuffled neighboring plants and by the genome rotation scheme (p > 0.1). Together with the significant PVE by neighbor genotype effects, these results suggest a heritable but polygenic basis for neighbor effects on the herbivore damage and insect community composition.

In the Neighbor GWAS, the sign of the estimated neighbor genotype effects represents positive or negative interactions between the two alleles of paired neighbors, which corresponds to associational resistance or susceptibility to herbivore damage and insect community composition. We therefore examined the number and frequency of SNPs with positive or negative separately for the four phenotypes per site. We detected both negative and positive in the top 0.1%-associated SNPs. The SNPs with associational resistance involved more minor alleles than those with associational susceptibility in five out of the four phenotypes × two sites, including the herbivore damage in Zurich. These data suggest that associational resistance may be less frequent than associational susceptibility among the A. thaliana genotypes. This frequency bias further motivated us to analyze another possible difference between the SNPs with positive and negative, such as patterns of selection scanned by using genome-wide polymorphism data. This selection scan revealed that SNPs with associational resistance had more signatures of balancing selection than those with associational susceptibility. These results highlight the different features of the top-scoring SNPs with associational resistance and susceptibility. Among many possible SNPs, we then attempted to determine key predictors that explain herbivore damage or insect community composition.

Genomic prediction and validation of key genotype pairs in the field
Due to the lack of significant SNPs associated with the neighbor effects, it was not feasible to predict key genotype pairs based on a few SNP predictors. We solved this problem using a genomic prediction approach that incorporated genome-wide SNPs for phenotype prediction. To predict the neighbor effects, we included all 1.2 million SNPs representing focal genotypes and neighbor genotypes in the least absolute shrinking and selection operator (LASSO). With or without neighbor genotypes, LASSO prediction was validated using a test dataset collected in an additional year. Among the four phenotypes we had measured per site, the test dataset of herbivore damage in Zurich was slightly better explained by the neighbor-including LASSO than by the neighbor-excluding LASSO (Spearman’s ρ = 0.416 and 0.391 for neighbor-including and neighbor-excluding LASSO, respectively). This result indicates that the herbivore damage in Zurich can be better predicted by incorporating neighbor genotypes. Therefore, we extracted 756 neighbor-related SNPs using the neighbor-including LASSO that better predicted the herbivore damage in Zurich.

To understand the potential mechanisms behind the 756 SNP predictors of the herbivore damage in Zurich, we performed gene ontology enrichment analyses for candidate genes relevant to associational resistance (the LASSO-selected SNPs with > 0) or associational susceptibility (< 0). For the candidate genes of associational resistance, we detected a significant enrichment of gene ontology with ‘jasmonic acid biosynthetic process’ (false discovery rate < 0.05), including the LIPOXYGENASE2 (LOX2) and LOX6 genes. LOX2 is particularly known as an essential gene for the production of volatiles, which can reduce herbivory on neighboring plants. In contrast, jasmonate-related gene ontology was not enriched in the candidate genes relevant to associational susceptibility. These results suggest the potential relevance of jasmonate-mediated defense in associational resistance to herbivore damage, motivating us to predict beneficial pairs using the set of SNPs selected by LASSO.

Using neighbor genotypes as a better predictor, we attempted to predict associational resistance or susceptibility to herbivore damage in Zurich. To this end, we extrapolated the 756 LASSO-selected SNPs to virtual mixture (a pair of two different genotypes) or monoculture (a pair of the same genotypes) conditions in silico. We estimated the relative effects of two-genotype mixtures on herbivore damage. Consistent with the higher frequencies of SNPs with associational susceptibility than associational resistance in GWAS, the pairwise effect size had a negative mode of distribution. Susceptible plant genotypes under monoculture imposed more damage on their counterparts when planted with another genotype, which was suggested by the negative correlation between the pairwise effect size and estimated herbivore damage under monoculture (Spearman’s ρ = –0.38; test for no correlation, p = 3.4E-08). In addition, based on the pairwise effect size of mixed planting, our simulations confirmed that the herbivore damage increased when the number of genotypes increased. These results suggest the prevalence of associational susceptibility to herbivore damage in Zurich among the 199 genotypes. In this situation, we asked whether identifying beneficial pairs would be feasible.

Despite the prevalence of negatively interacting pairs, 823 pairs had a positive estimate for mixed planting. To verify associational resistance at the stand level in situ, we planted three genotype pairs under monoculture and mixture conditions at the Zurich site. We tested Bg-2 and Uod-1 as a pair with a large positive effect; Vastervik and Jm-0 as a pair with a moderate positive effect; and Bro1-6 and Bla-1 as a pair with a slightly positive effect from the range of positive effect sizes. Consistent with this order of effect size, the pair of Bg-2 and Uod-1 indeed showed a significant reduction (average of 24.8% within a pair) in herbivore damage in the mixtures in the field (marginal means of the herbivore damage, t > 2.02, d.f. = 125, adjusted p < 0.05). The pair of Vastervik and Jm-0 also showed a significant reduction (average of 22.7% within a pair) in herbivore damage in the mixture compared with the average monocultures (marginal means of the herbivore damage, t > 2.22, d.f. = 125, adjusted p < 0.05). Expected from their smallest effect size, the pair of Bla-1 and Bro1-6 did not show a significant reduction in herbivore damage in the mixtures. Next, we conducted laboratory choice experiments in which black flea beetles were allowed to feed on the mixture of each of the three pairs. This experiment found significant differences in herbivore damage between Bg-2 and Uod-1 (likelihood ratio test, χ2 = 13.35, d.f. = 1, p < 0.001); and between Vastervik and Jm-0 (χ2 = 5.71, d.f. = 1, p < 0.05); but not between Bla-1 and Bro1-6 (χ2 = 0.87, d.f. = 1, p = 0.35), in agreement with observed effects of the field experiments. Taken together, field experiments and laboratory evidence demonstrated that candidate genotype pairs underpinned associational resistance to herbivory.

Methods
Field GWAS experiments
Plant genotypes. We used A. thaliana genotypes that were selfed and maintained as inbred lines, called “accessions.” To study the genomic variation responsible for biotic interactions, we overlapped our accessions with those used in the GWAS of microbial communities and glucosinolates. We used 199 accessions with a few additional accessions, all of which were genotyped in the RegMap and 1001 Genomes projects. Seeds of these accessions were obtained from the Arabidopsis Biological Resource Center (https://abrc.osu.edu/). The Santa-Clara accession was replaced with Fja1-1 in 2018 because the genotype of Santa-Clara was unavailable. For the genotype data, we downloaded a full imputed SNP matrix of 2029 accessions from the AraGWAS Catalog. Of the full 10,709,466 SNPs, we used 1,819,577 SNPs with minor allele frequency (MAF) > 0.05. Ler(gl1-1) and Col(gl1-2) were included as positive controls to test the known single-gene effects of GLABRA1 (GL1) on flea beetle resistance. The Ler or Col genome was assigned to the two gl1 mutants, with only the GL1 locus differing between the parental wild-type and gl1 mutants.

Field setting. To investigate the two distinct herbivore communities, we used field sites within or outside the natural distribution range of A. thaliana. As a native site, we used the outdoor garden of the University of Zurich-Irchel campus (Zurich, Switzerland: 47°23’N, 8°33’E, alt. ca. 500 m). As an exotic site, we used the Center for Ecological Research, Kyoto University (Otsu, Japan: 35°06’N, 134°56’E, alt. ca. 200 m). In the Otsu site, weeds were mown before the experiment, and the surroundings were covered with agricultural sheets before the experiment. In the Zurich site, each experimental block was placed in a separate bed that was not accessible to molluscan herbivores.

The field experiment at Otsu was conducted from late May to mid-June, and that at Zurich was conducted from late June to mid-July. The exact date of the field survey is annotated on the original data file. Plants were initially grown under controlled conditions and then planted in a field garden for three weeks. Seeds were sown on Jiffy-seven pots (33-mm diameter), and stratified under 4°C for a week. Seedlings were cultivated for 1.5 months under a short-day condition (8h light: 16h dark, 20°C). Plants were then separately potted in plastic pots (6 cm in diameter) filled with mixed soil of agricultural composts (Profi Substrat Classic CLED73, Einheitserde Co. in Zurich; Metro-mix 350, Sun Gro Co., USA in Otsu) and perlites at a 3:1 L ratio. The potted plants were covered with agricultural shading nets and acclimated to field conditions for a few days. A set of the 199 accessions and an additional Col-0 accession—namely, 200 individuals in total—was randomly assigned to each block without replacement and positioned in a checkered manner. Eight blocks of the 200 accessions were set at each site on 2017 and 2018 for GWAS. Three blocks of the 200 accessions were set at each site in 2019 for the model validation of LASSO (see “Modified Neighbor GWAS for LASSO” below). The blocks were > 2.0 m apart.

Phenotype survey. Insects and herbivorous collembola on individual plants were visually counted every 2–3 days. These species were identified using a magnifying glass. Dwelling traces and mummified aphids were also counted as proxies for the individual number of leaf miners and parasitoid wasps, respectively. Eggs, larvae, and adults were counted for all species, as long as they could be observed by the naked eye. All counts were performed by a single observer (Y. Sato) during the daytime at each site. At the Zurich site, yellow-striped and black flea beetles occurred every year (at least for four years). Small holes made by these flea beetles were counted at the Zurich site and their maximum number throughout the experiment was used as an indicator of herbivore damage. This phenotyping was not applicable to Otsu, because the most abundant herbivores were not flea beetles. Instead, the percentage of leaf area loss was scored in Otsu at the end of the experiment as follows: 0 for no visible damage, 1 for < 10%, 2 for > 10% and < 25%, 3 for > 25% and < 50%, 4 for > 50% and < 75%, and 5 for > 75% of area eaten.

We also recorded the initial plant size and presence/absence of inflorescences to incorporate these phenotypes as covariates in the statistical analyses. Initial plant size was evaluated by the length of the largest rosette leaf (mm) at the beginning of the field experiment because this parameter represents the plant size at the growth stage. The presence or absence of inflorescences was recorded 2 weeks after transplantation. Herbivore damage was evaluated by the number of leaf holes in Zurich and the leaf area loss in Otsu, as described above. The maximum number of individuals in each experiment was used as an index for the abundance of each insect species.

In our dataset, we defined extended plant phenotypes that represent insect community composition based on herbivore feeding habits and species richness. The insect community composition more significantly differed between the two sites than between 2017 and 2018 (redundancy analysis, F = 401, p < 0.001 for the sites; F = 152, p < 0.001 for the years); thus, the dataset was separated into Zurich and Otsu. The individual number of external chewers or other herbivores was defined as the total number of individuals of leaf-chewing species (e.g., beetles and caterpillars) or species that ate internal parts of a plant (e.g., phloem-sucking aphids, cell-content-sucking thrips, sap-sucking stink bugs, and leaf miners). The reason for this classification was to reflect the difference in plant defense responses to external chewers and other herbivores through jasmonic acid and salicylic acid pathways, respectively. Leaf miners are known as endophagous herbivore that can elicit both the jasmonic and salicylic acid pathways in plant defense responses. Based on their relevance to plant defense responses mediated by the salicylic acid pathway, leaf miners were included in the category of other herbivores. The individual number of other herbivores, to a large extent, represented the individual number of sap-sucking herbivores (i.e., phloem-sucking aphids, cell content-sucking thrips, and sap-sucking stink bugs) because leaf miners were very rare (0.35% of all individuals in the other herbivore category; see the raw data). Specialist–generalist classification was not applicable to our dataset because generalist herbivores were much fewer than specialist herbivores at both sites. Herbivore–carnivore ratio was also not applicable because carnivorous insects (e.g., parasitoid wasps and aphidophagous ladybirds) were much fewer than herbivores. These carnivorous insects were taken into consideration to calculate insect species diversity. For the index of insect species diversity, we calculated the exponential Shannon diversity and Simpson diversity indices in addition to the total number of insect species i.e., species richness. However, Shannon diversity and Simpson diversity showed a bimodal distribution that did not suit GWAS, and only the total number of insect species had quantitative phenotypic values. We, therefore, used the total number of insect species as an index of the insect species diversity. The insect community composition was analyzed using the vegan package v2.6-4 in R. All phenotypes except for the leaf area loss (score variable) were ln(x+1)-transformed to improve normality for GWAS and genomic prediction. Unless otherwise stated, all statistical tests were two-sided and all figure presentations and basic statistical analyses were performed using R version 3.6.1 or above.

GWAS with focal and neighbor genotype effects
Neighbor GWAS model. To separate the effects of focal and neighbor genotype effects on herbivore damage and insect community composition, we used the Neighbor GWAS model developed by Sato et al. (2021). Neighbor GWAS employs a linear mixed model that includes additional fixed and random effect. According to the same mixed model as Eq. 2 of Sato et al., we present a specific case that fits the data structure of the present field experiments as follows. Let xi and xj be the allelic status at each SNP of the i-th focal plant and j-th neighboring plant, respectively. The inbred accessions had two states as x ∈ {–1, +1}. A phenotype value of the i-th focal individual plant yi was then given as
yi = β0 + β1xi + β2(1/J)Σj=1J xixj + ui + ei (1)
where β0 is the intercept; β1xi is a fixed effect of the focal genotype and the same as standard GWAS, and the second coefficient β2 determines positive or negative effects from the mean allelic similarity (1/J)Σj=1J xixj at a given locus between the focal individual i and neighboring individuals j up to the total number of neighboring individuals J. When the focal plant shares the same allele with a neighboring plant, the product xixj = (–1) × (–1) = +1 or (+1) × (+1) = +1. By contrast, when a neighboring plant has a different allele from the focal plant, the product xixj = (–1) × (+1) = –1 or (+1) × (–1) = –1. The average of these products represent the mean allelic similarity of the focal plant with neighboring plants (1/J)Σj=1J xixj (J representing the number of neighboring plants L in the linear mixed model Eq. 2 of Sato et al.). The coefficient β2 then determines the direction and strength of the effects of neighbor genotype similarity on a phenotype. In addition to the fixed effects, the random effects ui and residuals ei contribute to phenotypic variation. These random effects and residuals follow a normal distribution as ui ∼ N(0, σ21K1 + σ22K2) and ei ∼ N(0, σ2e) (‘~’ meaning ‘distributed as’: see also page 3 in Sato et al.). The variance component parameters σ21 and σ22 represent the polygenic effects of focal and neighbor genotypes on a phenotype, respectively. Same as the standard GWAS, K1 represents a kinship matrix among n plants and is defined as the fraction of shared alleles among all SNP sites for a pair of accessions. K2 represents a sample structure due to additive effects of neighbor genotypes among n plants given by the cross-product K2 = XTX / (q–1) (see page 3, “Variation partitioning” in Sato et al.). The elements of X2 include covariates of the neighbor genotype similarity. With all the elemental covariates shown, the q-SNPs × n-plants matrix X2 can be expressed as follows:
X2 = [ (1/J)Σj=1J xq,1j, (1/J)Σj=1J xq,2j, ..., (1/J)Σj=1J xq,nj ]
where y is n × 1 vector of phenotypes; X includes mean, focal genotype values x, neighbor genotype similarity (1/J)Σj=1J xixj and other founding covariates as a matrix of fixed effects for n plants; ~β is a vector that represents coefficients of the fixed effects; Z is a design matrix that assigns individuals to a genotype, and becomes an identity matrix if all plants have different genotypes; u is the random effect with Var(u) = σ21K1 + σ22K2; and e is residual as Var(e) = σ2eI. We estimated σ^21 and σ^22 (‘hat’ meaning estimates) using Eq. 2, and calculated PVE by Neighbor GWAS model as PVE = (σ^21 + σ^22) / (σ^21 + σ^22 + σ^2e).
When β2 and σ22 are set to 0, the Neighbor GWAS model (Eq. 1) is equivalent to the standard GWAS model y = β0 + β1x + u + e or y = X~β + Zu + e, where Var(u) = σ21K1 (see page 4 of Sato et al.). We estimated σ^21 = σ^21 + σ^22 + σ^2e, where the standard GWAS model is equivalent to the self by Sato et al. Because focal genotype values xi and neighbor genotype similarity (1/J)Σj=1J xixj are correlated each other (Supplementary Fig. 14 and Supplementary Note 4), it was difficult to separate the focal and neighbor genotype effects as σ^21 = σ^21 + σ^22 + σ^2e and σ^22 = σ^21 + σ^22 + σ^2e (see also the page 4 and Fig. 3 in Sato et al.). Instead, the net contribution of neighbor genotype effects to a phenotype can be sufficiently evaluated as PVE – h2, which is equivalent to the net PVEnei of Sato et al. (Supplementary Fig. 15 and Supplementary Note 4). To avoid confounding focal and neighbor genotype effects, PVE should be tested from simpler to complex models. This stepwise procedure for the likelihood ratio tests was established in our previous study (see page 6 in Sato et al.), as follows:
1. Compute the null likelihood with σ21 = 0 and σ22 = 0.
2. Compare models with and without σ21 based on the likelihood ratio test for σ21.
3. Calculate h2 = σ^21 / (σ^21 + σ^22 + σ^2e).
4. Compute the likelihood with σ^21 ≠ 0 and σ22 = 0.
5. Compare models with and without σ22 based on the likelihood ratio test for σ22.
6. Calculate PVE = (σ^21 + σ^22) / (σ^21 + σ^22 + σ^2e).
This procedure was performed using the rNeighborGWAS package v1.2.4, where linear mixed models were solved using the average information-restricted maximum likelihood method implemented in the gaston package v1.5.5. The statistical significance of the PVE was determined using the likelihood ratio between a simpler and complex model, which asymptotically follows a χ2 distribution with one degree of freedom. Multiple testing was not considered among phenotypes as it was not common in PVE analyses and GWAS of multiple phenotypes. The initial plant size, presence/absence of inflorescences, study years, and differences in experimental blocks were considered non-genetic covariates. To empirically test the significance of PVE, we also performed a permutation test in which neighboring plants within blocks or genomic positions were randomized (Supplementary Fig. 6 and Supplementary Note 1). We repeated the GWAS analysis at J = 0 (i.e., standard GWAS), J = 4 (up to the nearest neighbors) and J = 12 (up to the second-nearest neighbors). If K2 is ignored, imperfect correction of the sample structure leads to inflation or deflation of the p-values (Supplementary Fig. 16 and Supplementary Note 5).

List of candidate genes. Provided that linkage disequilibrium (LD) decays within 10 kbp on average in the A. thaliana genome, we searched for candidate genes within 10 kbp near SNPs with the top 0.1% p-values. Functional annotation data from The Arabidopsis Information Resource (TAIR) were used for the gene model and description of A. thaliana.

LASSO with focal and neighbor genotype effects
Modified Neighbor GWAS for LASSO. To perform multiple regressions on all SNPs, we used sparse regression that could simultaneously select important SNP predictors and estimate their coefficients. The Neighbor GWAS model (Eq. 1) is expressed as a multiple regression model as follows:
y = X0β0 + X1β1 + X2β2 + e (3)
where y is a phenotype vector; β0 is a vector including coefficients for an intercept and non-genetic covariates; β1 and β2 are vectors including coefficients of focal and neighbor genotype effects, respectively; X0 is a matrix that includes a unit vector and non-genetic covariates for n individuals. X1 is a matrix that includes the focal genotype values for n individuals and q SNP markers. X2 is a matrix that includes the neighbor genotype similarity for n individuals and q SNP markers, as noted above (see the subsection ‘Neighbor GWAS model’ above). To simultaneously perform variable selection and coefficient estimation, we applied the least absolute shrinkage and selection operator (LASSO) to Eq. 3. We further cut off 1,819,577 GWAS SNPs to 1,242,128 SNPs for LASSO with the criterion of LD at r2 < 0.8 between adjacent SNPs, because LASSO is sensitive to high correlations among explanatory variables. The initial plant size, presence/absence of inflorescences, study years, and experimental blocks were considered fixed covariates. Important variables were selected from 1,242,128 SNP markers and the same number of neighbor-related SNPs using LASSO. We used the Python (v3.6.8) version of the glmnet package v1.0 to perform the LASSO. The kinship and sample structure among individuals were implicitly considered because LASSO regression can deal with all SNPs simultaneously. While a gradient of sparse regressions from the LASSO, via the elastic net, to the ridge regression was available in the glmnet package, we used the sparsest regression, LASSO, because of the computational burden of recursive calculation during the effect size estimation and simulation (see “The effect size of mixed planting” below).

To determine the LASSO regularization parameter λ, we first trained the LASSO models with the learning data (2017 and 2018) and then validated their outputs using the test dataset collected in another year (2019; see also “Field setting” above). The predictability of the four phenotypes was evaluated based on the correlations between the predicted and observed values of each phenotype. Spearman’s rank correlation ρ was used because some phenotypic values were not normally distributed. The predicted values were obtained from LASSO models with different values of λ. To assess genetically based predictability, we quantified the observed phenotype values in 2019 as the residuals of a standard linear model. This standard linear model incorporated the same non-genetic explanatory variables as the LASSO model, including the initial plant size, presence of inflorescences, and difference in three experimental blocks, while each phenotype was considered a response variable. To determine whether the incorporation of neighbor genotypes improved the correlation with the test data, we compared LASSO with or without neighbor genotypes across a series of λ values. If the neighbor-including LASSO yielded a larger correlation than the neighbor-excluding LASSO at a given λ, this indicates that neighbor genotypes were able to improve the predictability of a target phenotype by LASSO. In this context, the maximum ρ of the neighbor-including LASSO was larger than that of the neighbor-excluding LASSO for herbivore damage in Zurich (Supplementary Fig. 10a). Furthermore, the neighbor-including LASSO achieved this maximum ρ even with stringent regularization (= larger λ) compared to the neighbor-excluding LASSO (Supplementary Fig. 10a). For the Otsu site, the neighbor-including LASSO also had slightly larger correlations with herbivore damage than the neighbor-excluding LASSO, supporting the improved predictability of herbivore damage by neighbor genotypes at another site (Supplementary Fig. 10e). None of the community composition phenotypes, however, showed better predictability with neighbor-including LASSO (Supplementary Fig. 10b–d, f–h). This was presumably because the individual number of the predominant species differed between study years (Supplementary Fig. 1b–g). These additional results support the improved predictability of herbivore damage but suggest difficulty in predicting community composition by neighbor genotypes.

When the neighbor-including LASSO outperformed the neighbor-excluding ones at a given λ, we obtained the vectors of the estimated coefficients β^2 that were able to improve the phenotype prediction. LASSO could yield multiple sets of β^2 across a series of λ where the neighbor-including LASSO yielded larger correlations. A larger λ tends to yield fewer non-zero SNPs with large coefficients, whereas a smaller λ tends to yield more non-zero SNPs with small coefficients. To consider the polygenic basis of neighbor effects, we averaged the estimated coefficients β^2 per SNP across the range of λ, resulting in 756 SNPs with non-zero β^2 for herbivore damage in Zurich (see “Results” in the main text). This estimated vector of neighbor coefficients β^2 was used to estimate the effect size.

Post-LASSO analysis (i): The effect size of mixed planting. To estimate the pairwise effect size of mixed planting, we extrapolated the LASSO model (Eq. 3) under a virtual monoculture (a pair of the same accessions) or pairwise mixture (a pair of different accessions). The pairwise effect size was determined by the difference in the linear sum [xi°xj]β^2 – [xi°xi]β^2 between a pair of accessions. The first term [xi°xj]β^2 represents the phenotype values expected from different genotype vectors between accessions i and j (= pairwise mixture), whereas the second term [xi°xi]β^2 represents those expected from the same genotype vectors between accessions i and i (= monoculture). The element-wise product [xi°xj] or [xi°xi] represents the neighbor genotype similarity between a pair of different or the same accessions, respectively. The neighbor genotype effects turned out to have no significant SNPs (Figs. 2e–h, m–p, 3a–d, and i–l); therefore, the genotype pairs predicted by many moderate-effect loci were suitable for testing the estimated effects of mixed planting. In contrast, genotype pairs showing the largest effect size were selected based on a few large-effect but less reliable loci. Assuming that multiple moderate-effect loci could result in the effects of mixed planting, we avoided the extreme tail of the effect size distribution when focusing on the three pairs: a pair with a large positive effect, Bg-2 and Uod-1 (effect size = 0.8); a pair with a moderate positive effect, Vastervik and Jm-0 (0.23); a pair with a slight positive effect, Bro1-6 and Bla-1 (0.1) (Fig. 4a). Note also that β2 in the neighbor GWAS models (Eqs. 1 and 3) denotes symmetric interactions between the focal i and neighbor j individuals, and thus [xi°xj] and [xj°xi] exert the same effects on a target phenotype.

To test whether the increasing number of plant genotypes increases or decreases herbivore damage, we also simulated herbivore damage in Zurich, i.e., ln(no. of leaf holes + 1) using the estimated vector of neighbor coefficients β^2. Assuming the nearest neighbors in a two-dimensional lattice, we simulated mixtures of up to eight genotypes. The herbivore damage was predicated by its marginal value with respect to the net neighbor effects [xi°xj]β^2. To examine the overall and selected patterns, we tested two types of genotype selection: (i) random selection from all pairs or (ii) random selection from pairs with positive estimates of pairwise mixed planting (positive values in Fig. 4a). First, eight genotypes were randomly selected from the 199 accessions to represent the overall pattern (Fig. 4c). We listed one (monoculture), two, four, or eight (full mixture) genotype combinations among the selected eight genotypes and averaged their predicted damage [xi°xj]β^2 among all the combinations. Second, four positively interacting pairs (x-axis > 0 in Fig. 4a) were randomly selected to test whether the random selection of positive pairwise interactions yielded positive relationships between genotype number and anti-herbivore resistance (Supplementary Fig. 11f). Duplicates of accessions were not allowed when selecting four pairs of two paired accessions. This line of random sampling was performed 9999 times to calculate the mean and standard deviation. In the first case, Fig. 4c shows the negative relationship between the number of genotypes and plant resistance. In the second case, herbivore damage decreased by paired mixing but increased by four- and eight-genotype mixing (Supplementary Fig. 11f). This was because scaling up pairwise mixtures to four or eight genotypes confounded negatively interacting pairs. In addition to Fig. 4a, c, these supplementary results also support the difficulty in targeting the positive relationships between genotype richness and anti-herbivore resistance.

Post-LASSO analysis (ii): GO enrichment analysis. To infer the category of genes related to positive and negative neighbor effects, we performed gene ontology (GO) enrichment analyses for candidate genes near LASSO-selected SNPs (i.e., SNPs with non-zero β^2). Same as the list of candidate genes in GWAS, we searched for genes within 10 kbp around each selected SNP. We omitted duplicated genes after listing the candidate genes. We then performed Fisher’s exact probability tests for each GO category against the entire gene set of A. thaliana. Multiple testing for genes was corrected using a false discovery rate (FDR). The entire set was built upon TAIR GO slim annotation using the GO.db package v3.17.0 in R. To summarize the results of the GO enrichment analysis, we applied the REVIGO algorithm to the list of significant GO terms at FDR < 0:05. When summarizing the significant GO terms, we focused on the Biological Process with the similarity measure at 0.7 (i.e., the same as the default setting). The rrvgo v1.12.0 and org.At.tair.db v3.17.0 packages in R were used to run the REVIGO algorithm. This line of GO analysis was separately performed for SNPs with negative or positive β^2 to detect GO terms unique to positive or negative neighbor effects on anti-herbivore resistance. Note also that post-GWAS GO analyses possess the issue of statistical non-independence due to LD in the standard GWAS. However, LASSO was less likely to be subject to this issue because (i) this sparse regression could sparsely select SNP variables across a genome, (ii) we pruned adjacent SNPs on the strong LD at r2 > 0.8, and (iii) we focused on unique genes before using Fisher’s test. Therefore, we applied the conventional GO enrichment test based on Fisher’s test with FDR correction to the LASSO results. The in-house R package that includes utility functions of the GO enrichment analysis is available at Zenodo.

Mixed planting experiment
Field experiment. To test the effects of mixed planting on herbivore damage, we transplanted three pairs of accessions (i.e., Bg-2 and Uod-1; Vastervik and Jm-0; and Bla-1 and Bro1-6) under mixture and monoculture conditions. The theory of plant neighbor effects suggests that both the plant patch size and neighbor composition should be manipulated to distinguish the effects of mixed planting from the density-dependent attraction of herbivores. Therefore, we set large and small plant patches in addition to monoculture or mixture conditions. The field experiment was conducted from late June to July 2019 and 2021 in the outdoor garden of the University of Zurich-Irchel. Plants were first grown under short-day conditions and then transferred to the outdoor garden following the same procedure as the field experiment for GWAS above. Two accessions were then mixed in a checkered manner under the mixture condition, whereas either of the two accessions was placed under the monoculture condition. The large patch included 64 potted plants in 8 × 8 trays and had a single replicate, whereas the small patch included 16 plants in 4 × 4 trays and had three replicates (upper photographs in Supplementary Fig. 12). In the mixture setting, the two potted accessions filled the square space in a checkered manner without a blank position (upper photographs in Supplementary Fig. 12). The total number of initial plants was two accessions × three pairs × mixture or monoculture × large or small patches × two years = 2016 individuals. Only a few pots per plot were labeled to track the plots in the field, whereas the other pots were not labeled to blind their information. The initial plant size was measured in the same manner as in the field GWAS. Leaf holes were counted three weeks after transplantation. Four plants died during the field experiment, resulting in a final sample size of 2012 plants.

Statistical analysis. We analyzed the herbivore damage (i.e., the number of leaf holes per plant) as a response variable. Linear mixed models were used for the number of leaf holes because this variable appeared to be normally distributed. The number of leaf holes was ln(x+1)-transformed to improve the normality. The explanatory variables were plant accession, mixture or monoculture conditions, small or large patches, and study years. The initial plant size, represented by the length of the largest leaf (mm), was considered as an offset term. Two-way interactions were also considered among the plant accessions, mixture conditions, and patch conditions. Because the large and small patches had different numbers of individual plants, this imbalance was dealt with using a random factor. We split the large patch by 4 × 4 potted plants (= the same size as the small patch; see also photographs in Supplementary Fig. 12), and considered these subplot differences — i.e., the total of 126 subplots — as a random effect. The significance of each explanatory variable was tested using Type III analysis of variance based on Satterthwaite’s effective degrees of freedom and F-tests. To compare herbivore damage for each accession between the mixture and monoculture conditions, we calculated marginal means for the full model based on Satterthwaite’s method with Sidak correction of multiple testing for accessions. For the analyses of leaf holes, we used the lme4 v1.1-34, lmerTest v3.1-3, and emmeans v1.8.7 packages in R. Boxplots visualize the median with upper and lower quartile, with whiskers extending to 1.5 × inter-quartile range.
To examine the effects of patch size and year in addition to mixed planting (Fig. 4d), we analyzed a separate dataset for patch conditions and study years (Supplementary Fig. 12a–d and Supplementary Table 5). Consistent with the order of the estimated effect size (Fig. 4a), the marginal means across these conditions showed the largest sum of the effects of mixed planting between Bg-2 and Uod-1 (= 0.495 in Supplementary Table 4b) and the second largest effect between Vastervik and Jm-0 (= 0.453 in Supplementary Table 4b). The significant effects of mixed planting on herbivore damage were more detectable in the large patches than in the small patches (Supplementary Fig. 12). The Bg-2 and Uod-1 accessions showed a significant reduction in herbivore damage among five cases out of the two accessions × two years × two patch conditions (Supplementary Fig. 12 and Supplementary Table 5) and a marginally significant case in the small patch (p = 0.053 in Supplementary Table 5a). The Vastervik and Jm-0 showed three significantly positive cases favoring the reduction in herbivore damage out of the eight conditions (Supplementary Fig. 12 and Supplementary Table 5), indicating less consistency than the Bg-2 and Uod-1 pairs under diverse conditions. The Bla-1 and Bro1-6 pairs did not have significantly positive cases favoring the reduction in herbivore damage out of the eight conditions and even had one case of increased damage by mixed planting (Supplementary Fig. 12 and Supplementary Table 5). The main results and separate data show that the order of the observed mixing effects is consistent with that of the estimated effect size.

Laboratory choice experiment
Insect materials. To examine the feeding by flea beetles, we conducted laboratory choice experiments using one of the two major flea beetles, the black flea beetle Phyllotreta astrachanica. Adult P. astrachanica were collected from Brassica spp. at the University of Zurich-Irchel. Adults and larvae were reared on German turnips (Kohlrabi) following a previously established protocol. The species of flea beetles were identified based on the DNA sequence of the mitochondrial gene encoding cytochrome c oxidase subunit I (COI). DNA was extracted using ZYMO RESEARCH Quick-DNA Tissue/Insect Kits (cat. no. D6016). We used universal COI primers designed by Folmer et al. for Polymerase Chain Reaction (PCR) amplification under the following conditions: Initial denaturation at 95°C for 5 minutes followed by 40 cycles of 95°C for 15s, 50°C for 30s, 72°C for 60s and a final extension at 72°C for 3 min. The PCR products were sequenced using Sanger sequencing. We compared our sequences with the COI sequences registered by Hendrich et al., which included 15 Phyllotreta species with several individual vouchers per species collected in Central Europe. Our sequences and the registered sequences were clustered using a neighbor-joining tree and the default alignment method implemented in the Qiagen CLC Main Workbench. We identified the species from our samples based on phylogenetic clusters. Our sequence data were registered in GenBank with IDs OQ857829 to OQ857834, which included three individuals of black- and yellow-striped flea beetles.

Experimental setting. We used the three pairs of six A. thaliana accessions, Bg-2 vs. Uod-1, Vastervik vs. Jm-0, and Bla-1 vs. Bro1-6. Seeds were sown on Jiffy-seven pots (33-mm diameter) and stratified at 4 °C for a week. Seedlings were cultivated under long-day conditions (16h light: 8h dark, 22/20 °C) for 3 weeks, with liquid fertilizer added a week after the start of cultivation. Two adult beetles were allowed to feed on the mixture of two individuals × two accessions for three days under long-day conditions (Supplementary Fig. 13a). The feeding arena was constructed using a transparent plastic cup (129 mm in diameter and 60 mm in height) that enclosed four Jiffy-potted seedlings. Excluding cups without any infestation by P. astrachanica, we obtained 15–20 replicates of feeding arena per pair.

Statistical analysis. We analyzed the herbivore damage (i.e., the number of leaf holes per plant) as a response variable using generalized linear models. Negative binomial errors and a log-link function were chosen because the number of leaf holes was zero-truncated (Supplementary Fig. 13b). Plant accessions and arena IDs were included as the explanatory variables. Likelihood ratio tests based on a χ2-distribution were used after checking whether the ratio of residual deviance to the residual degree of freedom was nearly one. The significance of each explanatory variable was tested by excluding one variable from the full model. The glm.nb function in the MASS package in R was used for generalized linear models with negative binomial errors. Likelihood ratio tests showed that flea beetles showed a significant preference between Bg-2 and Uod-1 and between Vastervik and Jm-0 but not between Bla-1 and Bro1-6 (Supplementary Table 6). The effect of the experimental area on leaf holes explained deviance but was only significant in the Bg-2 and Uod-1 pairs (Supplementary Table 6).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Study by Sato et al. describes a novel approach to identifying plant genotype pairs which can confer associational resistance to herbivores. Using genome-wide association method the authors predicted candidate genotype pairs that could reduce herbivory through associational resistance. They then have been able to validate their predictions by planting several of such pairs and confirming associational resistance. Taken together, this study provides an interesting and promising proof of concept enabling to predict impacts of intraspecific mixed planting on herbivory and there is a potential to use such strategy in the future for crops. I do however feel that the authors take very narrow view of possible mechanisms of associational resistance and their study ignores many other mechanisms which can lead to associational resistance or associational susceptibility. At the very least, the existence of these other mechanisms need to be acknowledged, even if they were not included in the current study.\n\nMain comments:\n- The authors refer to associational resistance and associational susceptibility as underlying mechanisms of effects of mixed planting on herbivores (lines 50-51). This is incorrect, terms AR and AS are used to describe the phenomenon of reduced/increased herbivory in mixed planting, not the mechanisms, which can include e.g. reduced apparency of host plants in mixed stands, changes in quality of host plants for herbivores depending on their neighbours, higher predation and parasitism on herbivores in mixed stands etc.\n- It is also incorrect to claim that associational resistance and susceptibility are driven by plant-plant interactions among neighbouring individuals (lines 52-53). Plant-plant interactions between neighbours definitely play a role, but so do plant-herbivore and plant-natural enemy interactions. For instance, in the push-pull strategy of pest management associational resistance is achieved though manipulation of pest behaviour by planting neighbours which are repellent (push) or more attractive (pull) to herbivores than the crop plant, so through plant-herbivore and not plant-plant interactions.\n- Lines 53-57. Repellence of herbivores by neighbouring plants is only one of many possible mechanisms of associational resistance, as I mentioned above. Equally, ‘spill over’ effect is only one of many potential mechanisms behind associational susceptibility (e.g. some generalist herbivores feed on several host plants to minimize toxicity etc).\n- I found the authors' classification of herbivores into internal and external feeders is confusing. Normally herbivores such as leaf gallers, leaf miners and borers which feed inside plant tissues are called internal feeders, but the authors seem to use this term to refer largely to sucking insects (aphids and thrips) which feed externally. Classifying herbivores into chewing and sucking would make more sense as these herbivores are also known to trigger different plant defensive pathways (JA vs SA). There is also inconsistency in the text as to which herbivores are considered internal. In the methods (lines 412-413) it is stated that internal herbivores included aphids, thrips and leaf miners whereas in the main text (lines 83-84) only aphids and thrips are listed as internal herbivores.\n- Lines 112-113. Caterpillars are referred to as sedentary as they do not move between plants. However, plant selection is usually not done by the caterpillars but by the female moths/butterflies which lay eggs on the plant on which the caterpillars feed, and moths and butterflies are not sedentary and very much able to move (fly) between plants. I am therefore not sure whether the authors conclusion on lines 119-120 that “These patterns of herbivore damage and communities suggest that neighbor genotypes are more likely to influence mobile herbivores than sedentary herbivores” holds.\n- I am not sure what is meant by ‘the space of neighbor effects’ (line 475) – do you mean the spatial scale at which neighbor effect occurs?\n\nMinor comments:\n- Lines 369-371: 2019. “The field experiment at Otsu was conducted from late May to mid-June, and that at Otsu was conducted from late June to mid-July” – I assume it is a mistake and one of the experiments should be the Zurich one.", "reviewer_2": "The authors have analyzed the genetic effects of neighbors on herbivory damage in A. thaliana using GWAS. Mapping loci that interact with neighbors is a challenging and intriguing topic. If the model can function effectively, its potential applications are extensive. However, their neighbor GWAS seems statistically weak, and the demonstration of performance provided is insufficient. I would suggest that the authors consider improving the model by collaborating with experts in statistics and population genetics."}, "author_response": "Response to Reviewers\nItalic: original comments; Plain: our responses; Blue: manuscript texts\nTo Reviewer #1\nComment #1.\nStudy by Sato et al. describes a novel approach to identifying plant genotype pairs\nwhich can confer associational resistance to herbivores. Using genome-wide\nassociation method the authors predicted candidate genotype pairs that could\nreduce herbivory through associational resistance. They then have been able to\nvalidate their predictions by planting several of such pairs and confirming\nassociational resistance. Taken together, this study provides an interesting and\npromising proof of concept enabling to predict impacts of intraspecific mixed planting\non herbivory and there is a potential to use such strategy in the future for crops. I do\nhowever feel that the authors take very narrow view of possible mechanisms of\nassociational resistance and their study ignores many other mechanisms which can\nlead to associational resistance or associational susceptibility. At the very least, the\nexistence of these other mechanisms need to be acknowledged, even if they were\nnot included in the current study.\nResponse: We appreciate reviewer 1’s insightful comments on associational\nresistance (AR) and associational susceptibility (AS). We fully agree that AR and AS\nare not mechanisms but phenomena. The goal of our study was to predict such\n‘phenomena’ of AR or AS based on genome-wide polymorphism data. While\nmentioning other mechanisms in the revised Introduction and Discussions (see our\nreplies to the comments #2-6 below), we have more carefully discussed and made\nour goal clearer (Discussion, Line 289-315). Based on this standpoint, we would like\nto explain how the revised Introduction and Discussions concern AR and AS in\nresponse to the comments #2-6 below.\nComment #2.\nMain comments:\n- The authors refer to associational resistance and associational susceptibility as\nunderlying mechanisms of effects of mixed planting on herbivores (lines 50-51). This\nis incorrect, terms AR and AS are used to describe the phenomenon of\nreduced/increased herbivory in mixed planting, not the mechanisms, which can\ninclude e.g. reduced apparency of host plants in mixed stands, changes in quality of\nhost plants for herbivores depending on their neighbours, higher predation and\nparasitism on herbivores in mixed stands etc.\nResponse: We agree that AR and AS are not mechanisms and have now rewritten\nthem as ‘phenomena’ and introduced more examples including plant–herbivore,\nplant–plant, and plant–enemy interactions that may lead to AR and AS. The second\nparagraph of the revised introduction now reads “In anti-herbivore defense, both\npositive and negative effects of mixed planting on stand-level resistance have been\nreported in the literature7,11–13. These phenomena are referred to as associational\nresistance and associational susceptibility, respectively11,14,15. Associational\nresistance and susceptibility involve ecological interactions that are mutually non-\nexclusive, such as plant–herbivore, plant–plant, and plant–carnivore interactions11.\nIn plant–herbivore interactions, associational resistance or susceptibility occur when\n1\nchemical and physical plant traits jointly repel or attract herbivores to neighboring\nplants14–17. For instance, plant odor and apparency may affect the settlement of\nherbivores17–19, while physical barriers and toxic metabolites may alter herbivore\nbehavior and growth after settlement20,21. Plant–plant interactions may modulate the\nexpression of these plant traits through volatile-mediated communications22 or direct\ncompetition23. Plant–carnivore interactions may also lead to associational resistance\nor susceptibility when a mixture of plants compared with the component\nmonocultures attracts more or fewer natural enemies of herbivores24,25. The\ncomplexity of associational resistance and susceptibility makes it difficult to\ndistinguish between positively and negatively interacting genotype pairs for anti-\nherbivore resistance. Nonetheless, the application of associational resistance to\nagriculture is becoming more and more important to reduce the use of\ninsecticides17,26.” (Introduction, Line 48-64). New references include the following\ndiverse examples.\n[14] Tahvanainen, J. O. & Root, R. B. The influence of vegetational diversity on the\npopulation ecology of a specialized herbivore, Phyllotreta cruciferae (Coleoptera:\nChrysomelidae). Oecologia 10, 321–346 (1972).\n[15] White, J. A. & Whitham, T. G. Associational susceptibility of cottonwood to a box\nelder herbivore. Ecology 81, 1795–1803 (2000).\n[17] Pickett, J. A., Woodcock, C. M., Midega, C. A. & Khan, Z. R. Push–pull farming\nsystems. Current Opinion in Biotechnology 26, 125–132 (2014).\n[18] Hambäck, P. A., Agren, J. & Ericson, L. Associational resistance: Insect damage\nto purple loosestrife reduced in thickets of sweet gale. Ecology 81, 1784–1794\n(2000).\n[19] Jactel, H., Birgersson, G., Andersson, S. & Schlyter, F. Non-host volatiles\nmediate associational resistance to the pine processionary moth. Oecologia 166,\n703–711 (2011).\n[21] Bustos-Segura, C., Poelman, E. H., Reichelt, M., Gershenzon, J. & Gols, R.\nIntraspecific chemical diversity among neighbouring plants correlates positively with\nplant size and herbivore load but negatively with herbivore damage. Ecology Letters\n20, 87–97 (2017).\n[22] Karban, R. & Maron, J. The fitness consequences of interspecific eavesdropping\nbetween plants. Ecology 83, 1209–1213 (2002).\n[23] Barton, K. E. & Bowers, M. D. Neighbor species differentially alter resistance\nphenotypes in Plantago. Oecologia 150, 442–452 (2006).\n[22] Bernays, E. A. & Minkenberg, O. P. J. M. Insect herbivores: Different reasons for\nbeing a generalist. Ecology 78, 1157–1169 (1997).\n[24] Letourneau, D. K. The enemies hypothesis: Tritrophic interactions and\nvegetational diversity in tropical agroecosystems. Ecology 68, 1616–1622 (1987).\n[25] Aartsma, Y. et al. Spatial scale, neighbouring plants and variation in plant\nvolatiles interactively determine the strength of host–parasitoid relationships. Oikos\n129, 1429–1439 (2020).\n[26] Tooker, J. F. & Frank, S. D. Genotypically diverse cultivar mixtures for insect\npest management and increased crop yields. Journal of Applied Ecology 49, 974–\n985 (2012).\nComment #3.\n- It is also incorrect to claim that associational resistance and susceptibility are driven\nby plant-plant interactions among neighbouring individuals (lines 52-53). Plant-plant\n2\ninteractions between neighbours definitely play a role, but so do plant-herbivore and\nplant-natural enemy interactions. For instance, in the push-pull strategy of pest\nmanagement associational resistance is achieved though manipulation of pest\nbehaviour by planting neighbours which are repellent (push) or more attractive (pull)\nto herbivores than the crop plant, so through plant-herbivore and not plant-plant\ninteractions.\nResponse: We agree that not only plant–plant but also plant–herbivore and plant–\nenemy interactions play a crucial role in AR and AS. These three types of ecological\ninteractions are now introduced in the revised Introduction (see our reply to the\ncomment #2 above), where the push-pull system is cited as an example of plant-\nherbivore interactions (citation no. 17). According to the revised Introduction, we\ndiscuss our results about AR in more detail in the Discussion. In our study system,\nspecialist predators or parasitoids have not been observed for the major herbivores,\ni.e., Phyllotreta flea beetles; thus, we focus more on plant–herbivore and plant–plant\ninteractions than plant–enemy ones. The third paragraph of the revised Discussion\nnow reads “Associational resistance can occur through ecological interactions that\nare mutually non-exclusive, such as plant–herbivore, plant–plant, and plant–\ncarnivore interactions. Although elucidation of the detailed mechanisms is not the\ngoal of genomic prediction, our genome-wide analysis provides insights into the\npotential mechanisms of associational resistance. For instance, key SNPs near\nLOXs suggest a potential relationship between jasmonate-induced defense and\nassociational resistance to flea beetles, as LOX genes alter volatile production\nthrough jasmonic-acid pathways46,47. Such volatile chemicals are known to play three\nmultifunctional roles in associational resistance by repelling herbivores away from\nneighboring plants19, eliciting defense responses in neighboring plants22, and\nattracting carnivorous insects25. In our study system, carnivore-mediated interactions\nare unlikely to explain associational resistance because specialist predators or\nparasitoids have not been observed for flea beetles. This observation suggests that\nvolatile-mediated mechanisms allowing herbivore repellency or direct plant–plant\ncommunication could be two of the three ecological interactions in our study system.\nThe complex polygenic basis of neighbor effects remains, however, as a challenge in\nfurthering the mechanisms of associational resistance. Manipulating multiple genes51\nis ultimately needed to verify the set of small-effect genes responsible for\nquantitative variation in herbivore damage. To this end, our genome-wide study\npaves a way to identify candidate genes for further genetic studies.” (Discussion,\nLine 298-315).\nComment #4.\n- Lines 53-57. Repellence of herbivores by neighbouring plants is only one of many\npossible mechanisms of associational resistance, as I mentioned above. Equally,\n‘spill over’ effect is only one of many potential mechanisms behind associational\nsusceptibility (e.g. some generalist herbivores feed on several host plants to\nminimize toxicity etc).\nResponse: We have now incorporated more examples of AR and AS by classifying\nthem into three types of ecological interactions, as explained above (in response to\ncomments #2 and #3). We also recognize that diet mixing can drive AS but this\nbehavior is mainly reported for generalist herbivores (e.g., Bernays & Minkenberg\nagainst them, diet mixing likely was less relevant in our case. To keep the point\nclear, we would like to focus on the associational resistance to specialist herbivores\nin the revised manuscript as shown in our replies to the comments #2 and #3 above.\nComment #5.\n- I found the authors' classification of herbivores into internal and external feeders is\nconfusing. Normally herbivores such as leaf gallers, leaf miners and borers which\nfeed inside plant tissues are called internal feeders, but the authors seem to use this\nterm to refer largely to sucking insects (aphids and thrips) which feed externally.\nClassifying herbivores into chewing and sucking would make more sense as these\nherbivores are also known to trigger different plant defensive pathways (JA vs SA).\nThere is also inconsistency in the text as to which herbivores are considered internal.\nIn the methods (lines 412-413) it is stated that internal herbivores included aphids,\nthrips and leaf miners whereas in the main text (lines 83-84) only aphids and thrips\nare listed as internal herbivores.\nResponse: To clarify how and why the current classification was adopted, we have\nmoved text from the Methods to the Results. These details include three\nclarifications. First, the category of internal feeders included aphids, thrips, and leaf\nminers, but the leaf miners were rare. We have clarified this point in the Methods as\n“The individual number of external feeders to a large extent represented the\nindividual number of sap-sucking herbivores (i.e., phloem-sucking aphids, cell\ncontent-sucking thrips, and sap-sucking stink bugs) because leaf miners were rare\n(Supplementary Fig. 1; Supplementary Table 2).” (Methods, Line 463-466) and\nbriefly in the main Results as “…internal feeders (mostly aphids, thrips, and stink\nbugs but few leaf miners)…” (Results, Line 119-120). Second, we did not find any\nleaf gallers on A. thaliana in this study (see Supplementary Table 2). Third, we\nacknowledge in the initial Method section that chewers and suckers tend to trigger\nJA- and SA-mediated defense, respectively, and have now moved this point to the\nResults (Line 115-116). The whole text about these three points now reads as\nfollows “…Typically multiple insect species were found on each individual plant\n(average = 2.1 and 1.8 in Zurich and Otsu, respectively; Fig. 2d and l). The fact that\nthe major herbivore species observed in our study were specialists of Brassicaceae\n(Supplementary Fig. 1; Supplementary Table 2) led us to classify these insects into\nexternal and internal feeders, which are known to trigger different defense pathways\nmediated by jasmonic or salicylic acids39,40. To identify the insect functional groups\nresponsible for the herbivore damage, we quantified three extended phenotypes of\ninsect communities: the individual number of external feeders (flea beetles in Zurich\nor caterpillars in Otsu), the individual number of internal feeders (mostly aphids,\nthrips, and stink bugs but few leaf miners), and the total number of insect species per\nindividual plant (Fig. 2b-d and j-l; Supplementary Fig. 2)…” (Results, Line 111-121),\nwhere the citation numbers 39 and 40 refer to Mertens et al. (2021) Plant defence to\nsequential attack is adapted to prevalent herbivores. Nature Plants 7:1347-1353; and\nZüst, T. & Agrawal, A. A. Mechanisms and evolution of plant resistance to aphids.\nNature Plants 2, 15206 (2016), respectively.\nTo highlight insect communities, we moved one figure with histograms of\nphenotype values from the supplementary information to the main text (Figure 2).\nThe picture of a leaf mining trace is added to this Figure 2 to clarify that the category\nof internal feeders included leaf miners as well as sap-sucking herbivores.\n- Lines 112-113. Caterpillars are referred to as sedentary as they do not move\nbetween plants. However, plant selection is usually not done by the caterpillars but\nby the female moths/butterflies which lay eggs on the plant on which the caterpillars\nfeed, and moths and butterflies are not sedentary and very much able to move (fly)\nbetween plants. I am therefore not sure whether the authors conclusion on lines 119-\n120 that “These patterns of herbivore damage and communities suggest that\nneighbor genotypes are more likely to influence mobile herbivores than sedentary\nherbivores” holds.\nResponse: We added a new paragraph to discuss this point more carefully\n(Discussion, Line 317-326). We agree that the oviposition by mobile females results\nin herbivory by larvae but our experimental period, i.e. 3 weeks, seems short for\nmoths and butterflies to mature and disperse again on short-lived plants like A.\nthaliana. We discuss this point and what is expected in long-lived plant species, as\nfollows “…Based on detailed field surveys of the insect community, we found that\nsome herbivore species were influenced by neighboring plant genotypes while\nothers were not. Specifically, we detected a significant influence of neighbor\ngenotypes on jumping flea beetles but not on sedentary caterpillars. The initial\noccurrence of larval moths and butterflies could be determined through oviposition\nby mobile adult females to some extent, but long-lived plant species would be\nsuitable to observe the process in which these caterpillars become mature and\ndisperse again onto neighboring plants. When our method is applied to long-lived\nplants, it may enable us to predict impacts of intraspecific mixed planting on more\ndiverse herbivores at community level.” (Discussion, Line 317-326).\nComment #7.\n- I am not sure what is meant by ‘the space of neighbor effects’ (line 475) – do you\nmean the spatial scale at which neighbor effect occurs?\nResponse: Yes, we meant the spatial scale at which neighbor effect occurred. We\nhave reworded this point explicitly as “spatial scale” throughout the manuscript\n(Supplementary Note 4; Supplementary Figures 14 and 15).\nComment #8.\nMinor comments:\n- Lines 369-371: 2019. “The field experiment at Otsu was conducted from late May to\nmid-June, and that at Otsu was conducted from late June to mid-July” – I assume it\nis a mistake and one of the experiments should be the Zurich one.\nResponse: We thank the reviewer for pointing this out. We have corrected this point\nin the Methods (Line 418).\nTo Reviewer #2\nComment #9.\nA. thaliana using GWAS. Mapping loci that interact with neighbors is a challenging\nand intriguing topic. If the model can function effectively, its potential applications are\nextensive. However, their neighbor GWAS seems statistically weak, and the\ndemonstration of performance provided is insufficient. I would suggest that the\nauthors consider improving the model by collaborating with experts in statistics and\npopulation genetics.\nResponse: We appreciate the critical comments by reviewer 2 and have done the\nfollowing in response to the major concerns.\nWe expanded the description of methodological procedures in more detail\nthroughout the Introduction, Methods, and Supplementary Notes 1–3. These\nrevisions integrate the relevant points of our methods paper (Sato et al. 2021\nHeredity 126:597-614; see also the revised Methods), validate the incorporation of\ntwo kinship matrices (comment #10), and provide a permutation test (comment #11).\nWe believe that the revised manuscript can clarify the points of Neighbor GWAS.\nPlease see also the comments #10–14 below. The full methods and performance\ntests of Neighbor GWAS are described in our previous publications (Sato et al. 2021\nHeredity 126:597-614; Sato et al., 2023 Evolution 77:1145-57). The model of\nNeighbor GWAS employed in the present study is equivalent and a case-specific\nversion of Sato et al. (2021). To highlight this point, the specific equations, pages,\nand figures of the previous publications are mentioned in the revised Methods and\nSupplementary Notes when necessary.\nIn Sato et al. (2021) we had designed and validated the model of Neighbor\nGWAS with an expert of quantitative genetics, Dr. Eiji Yamamoto, who has\npublications in GWAS, plant breeding, and gene-by-environmental modeling (e.g.,\nYano & Yamamoto et al. 2016 Nature Genetics; Yamamoto et al. 2017 Heredity;\nYamamoto et al. 2021 G3). The collaboration with this expert enabled us to deal with\nthe issue of the underlying variance structure (see also our reply to comment #10). In\nSato et al. (2021), we also validated the model performance using simulated data.\nThese simulated data showed conditions under which the Neighbor GWAS model\neffectively functions. For instance, there we tested how the performance of Neighbor\nGWAS depends on spatial scale and demonstrated that when a reference spatial\nscale is narrow and matches the true effective spatial scale of neighbor effects,\ncausal SNPs are detected very well (see figure reproduced below with left panels\nshowing median AUC > 0.9 for the case of distance decay parameter, α = 3,\nPVEself:PVEnei = 1:8, and spatial scale s = 1). This case represents the short-range\nneighbor effects at J = 4 in the present study, suggesting that non-significant SNPs\nwere not due to the model but due to other factors such as a polygenic basis of our\ntarget traits.\nestablished through previous publications. To address specific concerns raised on\nthe current manuscript, we have partly adapted the models and methods of our\nprevious papers as explained in our replies to comments #10–14 below.\nComment #10.\nRegarding the model, I do not understand why the authors needed to use two\nkinship matrices. The neighbor effects (covariates synthesized by the neighbor’s\ngenotype) are mostly corrected by K_1. Also, if I understand correctly, neighbor\nGWAS might be based on specific hypotheses, such as the chemical or biomass\neffects from neighboring plants. While strong chemical effects may be mappable,\nbiomass effects cannot be mapped due to their polygenic nature (as shown in Fig\n1B). Have the authors considered using variance component models with K_1 and\nK_2, instead of genome partitioning way (Yang et al., 2011, Nat.\nGenet; https://www.nature.com/articles/ng.823)?\nResponse: We considered K_2 to separate SNP-marker effects from polygenic\nvariance due to covariates synthesized by the neighboring genotypes. Such\nadditional random effects should be incorporated not only for Neighbor GWAS but\nalso for other GWAS methods that include multiple fixed effects (e.g., GWAS of\nindirect genetic effects by Montazeaud et al. 2023 Nature Ecol & Evol). Our previous\nsimulation for Neighbor GWAS agreed with this notation and showed that the\nincorporation of K_2 yielded better performance to detect causal variants (mean\ndifference of the area under ROC curve = 0.014 between models with or without\nK_2, p-value<10-16 by paired t-test among all simulated phenotypes; see page 11 in\nSato et al. 2021). In the present study, we provide additional analyses and now show\nthat the imperfect correction of population structure leads to the inflation of p-values\nunless K_2 is considered. The incorporation of additional random effects is therefore\nnecessary to avoid false positive results. We show this point using simulated data\nand real data about the herbivore damage observed at the Zurich site\n(Supplementary Note 5; Supplementary Figure 16). The Supplementary Note 5\n7\ndescribes these results as follows “Supplementary Note 5. Underlying sample\nstructure in Neighbor GWAS\nIn the main text, we included both 𝐊 and 𝐊 to correct for the sample structure in\n! \"\nboth the focal and neighbor genotype effects on a phenotype. Here, we excluded 𝐊\n\"\nto test whether 𝐊 alone was sufficient for correcting the sample structure among\n!\nindividuals. We analyzed the simulated data from Supplementary Note 4 with or\nwithout 𝐊 . To calculate the slope of the QQ plots with or without 𝐊 , we estimated\n\" \"\nthe neighbor genotype effects 𝛽 at the true spatial scale at 𝐽 = 8 (corresponding to\n\"\nthe nearest neighbors) in the simulated space. Analysis of the simulated data\nshowed that the 𝑝-values of 𝛽 were inflated unless 𝐊 was incorporated as an\n\" \"\nadditional random effect (Supplementary Fig. 16a and b). The inflation of 𝑝-values\nwas attributable to the imperfect correction of the sample structure due to 𝐊 ,\n\"\nbecause the simulated data did not confound environmental factors. To test the\ninfluence of 𝐊 using real data, we also analyzed the herbivore damage in Zurich at 𝐽\n\"\n= 4 without 𝐊 . This additional analysis of read data showed that the 𝑝-values of 𝛽\n\" \"\nwere inflated when 𝐊 was not considered (Supplementary Fig. 16c). These results\n\"\nstrongly suggest that 𝐊 should be included to avoid false positive results.”\n\"\nSupplementary Figure 16 visualizes these results, as follows.\nSupplementary Figure 16. Inflation of 𝑝-values in the simulated data and real\ndata on herbivore damage when 𝐊 was not considered. a Box plots comparing\n\"\nthe slope coefficient of the quantile-quantile (QQ) plots when 30 simulated\nphenotypes were analyzed using the mixed models without 𝐊 . The deviation from\n\"\none (a horizontal dashed line) indicates inflation or deflation of 𝑝-values. b An\nexample of QQ plots comparing the results of the analysis of simulated data with\n(blue) or without (red) 𝐊 . c A QQ plot showing observed and expected 𝑝-values for\n\"\nthe neighbor genotype effect 𝛽 on herbivore damage (leaf holes) in Zurich at 𝐽 = 4\n\"\nwithout 𝐊 . The diagonal dashed line indicates that 𝑦 = 𝑥.\n\"\nThe way to generate simulation data followed that of our methods paper (Sato et al.\n2021) and is now described in Supplementary Note 4. Please see also our response\nto comment #12.\nRegarding the estimation of PVE in previous Figure 2, we did not perform\ngenome partitioning per chromosome (Yang et al. 2011) but variance component\nestimation with K_1 and K_2. The calculation of PVE is now detailed in the Method\nsection (see the subsections “Phenotypic variation explained by focal and neighbor\ngenotype effects” and “Genome-wide association study (GWAS)”). In addition to the\nMethod section, the spatial-scale dependency of these PVEs is shown in\nSupplementary Figure 15b and Supplementary Note 4 using the simulated data (see\nalso our response to comment #12 below).\nAdditionally, have the authors conducted permutation tests for Neighbor GWAS?\nI believe it's essential to empirically validate \\beta_2 by random sampling to establish\nits distribution.\nResponse: We thank the reviewer for this comment and have conducted\npermutation tests to empirically determine the genome-wide threshold for β .\npermutations, (i) using the first two blocks of the present data for 99 permutations\nand (ii) using the full dataset of leaf holes in Zurich for 20 permutations. In both\ncases, neighboring plants were randomly shuffled within blocks as suggested by\nreviewer 3 (see also our response to comment #18). The permutation test of the\nsubset and full data both did not find any SNPs above the empirical threshold for β ,\n2\nbut the observed PVE remained far more than 5% significant among the permuted\nPVEs (Supplementary Figure 6). These supplementary results are consistent with\nthose of the main text reporting the significant PVE without non-significant SNPs,\nnamely a polygenic basis. We have now incorporated these additional results in the\nSupplementary Note 1 and Supplementary Figure 6 as follows “Supplementary\nNote 1. Permutation test for the genome-wide significance threshold\nBonferroni correction is often considered conservative in GWAS because\nlinkage disequilibrium (LD) makes individual SNPs not statistically independent85,86.\nPermutation tests solve this issue by empirically determining the significant\nthresholds, although they are computationally demanding86. Using permutation tests,\nwe attempted to determine the significance level for the proportion of phenotypic\nvariation explained (PVE) by neighbor effects and GWAS. To test whether the\nobserved PVE was more significant than random, we randomly permuted\nneighboring plants within blocks and calculated the PVE for the herbivore damage in\nZurich at 𝐽 = 4. The net contribution of neighbor genotype effects was determined as\nPVE by the Neighbor GWAS model – the genomic heritability ℎ\" (see Methods\n“Phenotypic variation explained by focal and neighbor genotype effects”). We\ncalculated PVEs for 99 permutations and tested whether the observed PVE was\nhigher than the 95 percentile of the permuted PVEs. In contrast to PVE, the full data\nfor GWAS permutations required a computational time of approximately 13 h per\niteration per phenotype in a 24-thread parallel calculation. To deal with this difficulty,\nwe conducted a large and small set of permutations using the full data and subset\ndata of the herbivore damage in Zurich at the number of neighboring individuals 𝐽 =\n4. For the full data, we permuted neighboring plants within blocks 20 times and\nobtained the maximum –log (𝑝) among the 20 iterations. For the subset data, we\n10\nfocused on the data collected from the first two blocks (= 200 individuals × 2 blocks)\nin the Zurich site in 2017 because the highest SNP observed in these two blocks (–\nlog (𝑝) = 6.53) was very similar to that observed among all the phenotypes in this\n10\nstudy (–log (𝑝) = 6.55 for the internal feeders in Otsu at 𝐽 = 4; Fig. 3k;\nneighbor genotype effects 𝛽* was obtained per iteration. Their 95 percentile was\n\"\nconsidered an empirical significance level of 𝑝 = 0.05 for the subset permutation.\nAs a result of permutations, we found that the observed PVE by neighbor\ngenotype effects on the herbivore damage was higher than the 95 percentile of\npermuted PVEs in Zurich (Supplementary Fig. 6a), supporting the significance of\nPVE by neighbor genotype effects. When Neighbor GWAS at 𝐽 = 4 was applied to\npermuted full data, we found a higher association score of permuted –log (𝑝) than\n10\nthat of all tested SNPs for the herbivore damage in Zurich (max. –log (𝑝) = 7.60\n10\namong 20 permutations: Supplementary Fig. 6b). To further examine the genome-\nwide significance level of GWAS, we analyzed the subset data of herbivore damage\nobserved in Zurich with 99-times permutations of neighboring plants. This additional\npermutation on the subset data determined the 5% empirical threshold of –log (𝑝)\n10\nas ‘6.738’ at 𝐽 = 4 (Supplementary Fig. 6c). This empirical threshold is close to the\nBonferroni threshold at 𝑝 = 0.3, namely –log (𝑝) = ‘6.783’ in our case of 1,819,577\n10\nnone of the four phenotypes per site had more significant SNPs at 𝐽 = 4 than 6.783.\nThese supplementary results provide consistent evidence for significant PVE without\nSNPs above the genome-wide threshold, suggesting a polygenic basis of neighbor\ngenotype effects.” Supplementary Figure 6 visualizes these results as follows.\nSupplementary Figure 6. Permutation tests for the genome-wide significance\nlevel and the proportion of phenotypic variation explained (PVE) by neighbor\ngenotype effects. a Histogram of 99-times permuted PVEs for the full data on\nherbivore damage in Zurich at 𝐽 = 4. PVE by neighbor genotype effects was\ndetermined as PVE by the neighbor GWAS model – the genomic heritability ℎ\",\nwhich showed the net contribution of neighbor genotype effects to phenotypic\nvariation. The vertical solid and dashed lines indicate the observed PVE and the 95\npercentile of the permuted PVEs, respectively. b 20-times permutation tests of the\nneighbor genotypic effects 𝛽 using the full data on the herbivore damage (leaf\n\"\nholes) in Zurich at 𝐽 = 4. The horizontal black line indicates the Bonferroni threshold\nat 𝑝 = 0.05. The horizontal gray line indicates the maximum –log (𝑝) among the\n10\npermutations. c 99-times permutation tests of the neighbor genotypic effects 𝛽\n\"\nusing the subset data on herbivore damage in Zurich at 𝐽 = 4. The horizontal black\nline indicates the Bonferroni threshold at 𝑝 = 0.05. The horizontal gray line indicates\nthe permuted threshold at 𝑝 = 0.05.\nWe could not perform the permutation test for all the phenotypes because a\nsingle iteration of Neighbor GWAS requires approximately 13 hours per iteration per\nphenotype using muti-threaded computation. A possible way to speed up GWAS\nmay be to apply gBLUP-based single-shot GWAS (e.g., Gualdrón Duarte et al. 2014\nBMC Bioinformatics; Bernal Rubio et al. 2016 Animal Genetics; Chen et al. 2017\nGenetics) per iteration; however, it remains mathematically unknown whether this\ntechnique can be applied to the case of Neighbor GWAS that includes two fixed\neffects and two random effects. The application of this technique might be an outlook\nto make Neighbor GWAS more efficient but is beyond the scope of the present\nempirical work.\nComment #12.\nCould the authors explain why p-values inflate as J increases? (See Fig S4).\nlinearity between the focal and neighbor SNP variables. Figure S2 in Sato et al.\n(2021), reproduced below, exemplifies how the focal and neighbor SNP variables are\ncorrelated in response to changes in minor allele frequency (MAF) and changes in\nspatial scale of neighbor effects.\nBy expanding this previous work, we now show how such co-linearity leads to\ninflation/deflation of p-values and reach the conclusion that ‘the reference spatial\nscales should not be broadened beyond the spatial scale at which PVE becomes\nsignificant’ (Supplementary Note 4) to avoid the risk of false positives. To show this\npoint, we performed additional simulations and analyses in Supplementary Note 4 as\nfollows: “Supplementary Note 4. Spatial scale of Neighbor GWAS\nWe showed the results at the spatial scale that encompassed the nearest neighbors,\nnamely at the number of neighboring individuals 𝐽 = 4, in the main text. Focusing on\ncorrelated variables between focal and neighbor genotype effects, we address why\nthe reference spatial scales should not be broadened beyond the spatial scale at\nwhich PVE becomes significant. As described in the page 3 of Sato et al. (2021)34,\nNeighbor GWAS can suffer from the model structure where multiple explanatory\nvariables with the focal genotypic value 𝑥 and neighbor genotype similarity\n#\n,∑%\n𝑥 𝑥 ./𝐽 are correlated in a single regression model (sensu collinearity). This\n$&! # $\nissue matters not only for Neighbor GWAS but also for other models that include\ngene-by-gene components90,91. In such models, the covariance of the kinship matrix\n𝐊 , which represents the genetic similarity between individuals, helps fit the\n!\nphenotypic variation attributable to gene-by-gene interactions. In the Neighbor\nGWAS, the collinearity problem can become severe because of minor allele\nfrequency (MAF) and the spatial scale 𝐽 (Supplementary Fig. 14; see also Figure S2\nof Sato et al. 202134). This occurs because, when MAF is low, neighbors 𝑥 are less\n$\nlikely to vary in space and most plants have similar values for neighbor genotype\nsimilarity\n,∑%\n𝑥 𝑥 ./𝐽. Furthermore, if the neighbor effects range was too broad to\n# $&! # $\nencompass an entire field i.e., 𝐽 → ∞, the neighbor covariate and focal genotype\nvalue 𝑥 would become collinear as follows:\n,∑%\n𝑥 𝑥 ./𝐽 = 𝑥\n,∑%\n𝑥 ./𝐽 = 𝑥 𝑥‾ ,\n# $&! # $ # $&! $ # $\nwhere 𝑥‾ indicates a population-mean of neighbor genotypes and thus corresponds\n$\nto a population-mean of focal genotype values 𝑥‾ when 𝐽 → ∞. Two problematic\n#\noutcomes may occur when collinearity is severe. First the effect size of neighbor\ngenotypes (i.e., 𝛽 divided by the variance of 𝛽 estimates) may be overestimated if\n\" \"\nfocal genotypic effects explain most variations and thus underestimate the variance\nfor 𝛽 . Second, the neighbor genotype effects 𝛽 may be underestimated if\n\" \"\ncollinearity makes it unstable and overestimates the variance of 𝛽 . Below, we show\n\"\nhow these issues concern the simulated and real data.\n4.1. Simulation\nWe applied Neighbor GWAS for simulated data to examine the risk of false positive\nor negative signals with increasing 𝐽. Our previous publication showed that the\nNeighbor GWAS model accurately estimated PVE and effectively detected causal\nSNPs using extensive sets of simulation parameters (see the page 5 and Figure 4 in\nSato et al. 202134). In this study, we evaluated the inflation or deflation of 𝑝-values\nbased on quantile-quantile (QQ) plots. According to our previous study (3–5 in Sato\net al. 202134), we simulated polygenic traits in a rectangular space that assumed the\nsetting of the present field experiments using A. thaliana.\nGenotype: We used part of the A. thaliana RegMap panel32 to consider a realistic\ngenetic structure in the simulation. The genotype data on 1,307 accessions were\navailable on the Joy Bergelson laboratory website\n(http://bergelson.uchicago.edu/?page_id=790, accessed February 9, 2017). We\nextracted the first and second chromosome data with MAF >0.1, providing a SNP\nmatrix of 1307 plants with 65,226 SNPs. To simulate a rectangular space similar to\nthat in our field experiment using A. thaliana (Fig. 1), we randomly selected 324\naccessions (9 × 36 accessions) without any replacement and placed them in a 9 ×\n36 lattice space for each iteration.\nPhenotype simulation: We used Eq. 1 in Methods to generate phenotypes based on\nspatially arranged genotypes. We controlled the proportion of phenotypic variation\nwith the total PVE by genetic factors at 0.8, in which 0.6 out of 0.8 was attributed to\nSNPs with nonzero fixed effects, and residuals at 0.2. To simulate a polygenic\nphenotype, the number of causal SNPs was set to 50. The nonzero 𝛽 for the 50\n\"\ncausal SNPs was randomly sampled from −1 or 1 digit. The true number of\nneighbors was set at 𝐽 = 8 (up to the first nearest neighbors) to simulate the case of\nherbivore damage in Zurich in the main text. These parameter settings provide\nconditions under which the performance of the Neighbor GWAS model differentiates\nwell34. We used the nei_simu function implemented in the rNeighborGWAS\npackage34 to simulate the phenotypes. We ran 30 iterations of the genotype\narrangement and phenotype simulation.\nAnalysis of simulated data: We evaluated the inflation or deflation of 𝑝-values by\nfitting a linear regression to the QQ plots of the GWAS. The deviation of the slope\ncoefficient from one indicates inflation or deflation of the observed 𝑝-values than\nrandom expectation. The present simulation showed a larger deviation in the slope\ncoefficient when 𝐽 became larger than the true spatial scale, indicating the risk of\nfalse positive or negative signals (Supplementary Fig. 15a). In addition, we\ncalculated the net contribution of neighbor genotype effects, i.e., PVE – ℎ\", which is\nequivalent to the net PVE of Sato et al. (2021)34. This additional analysis showed\nnei\nthat the net PVE saturated at the true spatial scale at 𝐽 = 8 (Supplementary Fig.\nnei\n15b), which was consistent with the conclusion of our previous study34. These results\nfrom the PVE and QQ plot analyses indicate that the reference spatial scale should\nnot be broadened beyond the spatial scale at which the PVE is saturated.\n4.2. Long-range neighbor effects at 𝑱 = 12 in real data\n= 4 (main text) and 𝐽 = 12 (Supplementary Figs. 3-5) using real data. The reference\nspatial scale of the neighbor effects was broadened to the second-nearest neighbors\ni.e., 𝐽 = 12 for the herbivore damage and three extended phenotypes. Consistent\nwith the simulation, phenotypes that were significantly influenced by the nearest\nneighbors at 𝐽 = 4 exhibited an inflation of 𝑝-values at 𝐽 = 12 (Supplementary Figs. 4\nand 5), indicating a risk of false positives. For example, these inflated 𝑝-values at 𝐽 =\n12 were observed for the herbivore damage in Zurich and Otsu (Supplementary\nFigs. 4i and 5i); and for the individual number of external feeders in Zurich\n(Supplementary Figs. 4j). These results suggest that the reference spatial scales\nshould not be broadened beyond the spatial scale at which PVE becomes\nsignificant. To avoid false positive signals, we have presented the results of 𝐽 = 4 in\nthe main text, while including the results of 𝐽 = 12 for PVE (Supplementary Fig. 3),\nGWAS (Supplementary Figs. 4i-l and 5i-l), and other relevant analyses\n(Supplementary Figs. 8-10) in the Supplementary Information.” The data simulated\nduring this Supplementary Note 4 were also analyzed to address why K_2 should be\nincorporated in the mixed models (see the comment #10 above).\nComment #13.\nI couldn't find the motivation for the authors to begin testing the selection signature. It\nseems unnatural for A. thaliana to have different populations growing next to each\nother. Natural selection may not be necessary to influence 'neighbor A. thaliana.'\nResponse: As reviewer 3 pointed out, the past signature of selection we analyzed\nhere does not necessarily match what happens in the present field experiment.\nNatural accessions might not experience neighbor effects in an evolutionary time\nscale because they grew among different populations. Instead of this past selection\npattern, we focused more on the present patterns of β in our GWAS analysis. The\n2\nfurther examination of SNPs with positive or negative β found that associational\n2\nresistance may be less frequent than associational resistance in terms of minor allele\nfrequency (new Figure 3).\nIt is still possible to some extent that neighbor genotype effects occur within a\nnatural population of A. thaliana because within-population polymorphisms are\nreported for natural populations of A. thaliana and thus within-population GWAS is\nfeasible (e.g., Gloss et al. 2022 Phil Trans R Soc B). As reviewer 3 finds the\nselection scan useful (see comment #17), we would like to keep it in the\nSupplementary Information (Supplementary Note 4 and Supplementary Figure 11).\nBalancing between reviewer 2 and 3’s comments, the revised main text now\nfocuses more on β and puts less emphasis on the selection scan. The revised\ngenotype effects 𝛽* represents positive or negative interactions between the two\n\"\nalleles of paired neighbors, which corresponds to associational resistance or\nsusceptibility to herbivore damage and insect community composition34,43 (right inset in\nFig. 1a; Supplementary Fig. 7c-f). We therefore examined the number and\nfrequency of SNPs with positive or negative 𝛽* separately for the four phenotypes\n\"\nper site. We detected both negative and positive 𝛽* in the top 0.1%-associated SNPs\n\"\n(histograms in Fig. 3e-h and m-p; Supplementary Fig. 8a-h). The SNPs with\nassociational resistance (𝛽* > 0) involved more minor alleles than those with\n\"\nassociational susceptibility (𝛽* < 0) in five out of the four phenotypes × two sites\n\"\n13\n(boxplots in Fig. 3e-h and m-p; Supplementary Fig. 8i-p), including the herbivore\ndamage in Zurich (Fig. 3e). These data suggest that associational resistance may be\nless frequent among the A. thaliana genotypes. This frequency bias further\nmotivated us to analyze another possible difference between the SNPs with positive\nand negative 𝛽* , such as patterns of selection scanned by using genome-wide\n\"\npolymorphism data. This selection scan revealed that SNPs with associational\nresistance (𝛽* > 0) had more signatures of balancing selection than those with\n\"\nassociational susceptibility (𝛽* < 0) (Supplementary Fig. 9; Supplementary Note 2).\n\"\nThese results highlight the different features of the top-scoring SNPs with\nassociational resistance and susceptibility. Among many possible SNPs, we then\nattempted to determine key predictors that explain the herbivore damage or insect\ncommunity composition.” (Results, Line 187-205).\nComment #14.\nIn general, the statistical details and descriptions appear to be somewhat lacking.\nFor instance, on L92, please provide a thorough and precise description of the\n'significant SNP,' including its p-value, minor allele frequency (MAF), peak position,\nand its linkage disequilibrium (LD) information. What evidence supports the authors'\nclaim regarding the SNP's causal effect on the function of the GLABRA1 gene?\nResponse: Statistical details were given in the Methods or Supplementary Tables in\nthe initial version. We have now moved part of this information to the main text for\nfurther clarity. These descriptions are now partly redundant between main text and\nSupplementary Information, but we hope that a relaxed length limit now allows us to\naccommodate important information for the following four points.\nFirst, we moved necessary statistics from Supplementary Tables to the main\ntext like “likelihood ratio test, χ2 = xx, d.f.=1, p < 0.05”, “p < 0.05 after Bonferroni\ncorrection”, or “marginal means of herbivore damage, t = xx, d.f. = yy, adjusted p <\n0.05.” (Results, Lines 127, 132, 149, 153, 163, 171, 174, 182, 260, and 263).\nSecond, the significant SNP in GWAS was determined by the genome-wide\nBonferroni threshold at p<0.05 as previously described in the legend of Figure 3. We\nhave clarified this point again by adding the following sentence “The genome-wide\nsignificance level for 𝑝-values of 𝛽 or 𝛽 was determined using the Bonferroni\n! \"\ncorrection of multiple testing for all SNPs at 𝑝 = 0.05” to the Methods section (Line\n609-611).\nThird, as previously described in the Methods section (currently Line 405-\n408), we included two gl1 mutants as a positive control to test if our field GWAS\ndetected known effects of trichomes on herbivore damage. Our previous work at the\nsame field site reported the single-gene effects of the GL1 mutation on flea beetle\nresistance (Sato et al. 2019 BMC Biology). For clarity, we now shortly describe this\npoint in the Results as well as the Methods, as follows “Our previous study detected\nsingle-gene effects of the trichome developmental gene GLABRA1 (GL1) on\nresistance to herbivore damage made by flea beetles41; thus, two glabrous mutants\nwere included to test this effect. As expected, we detected significant effects of the\nmutation in the GL1 gene on the herbivore damage in Zurich (𝑝 < 0.05 after\nBonferroni correction; the most significant SNP being located on the third\nchromosome in Supplementary Fig. 4a; Supplementary Data 1).” (Results, Line 128-\n134).\nSNPs near LOX2 and LOX6 in the main text as follows “For the candidate genes of\nassociational resistance, we detected a significant enrichment of gene ontology with\n‘jasmonic acid biosynthetic process’ (false discovery rate < 0.05; Supplementary\nTable 4a), including the LIPOXIGENASE2 (LOX2) and LOX6 genes (chr3-16519704\nwith MAF = 0.40 near LOX2; chr1-25316686 with MAF = 0.44 and chr1-25317359\nwith MAF = 0.44 near LOX6; Supplementary Data 2).” (Results, Line 224-229). As\npreviously described in the subsection of the Methods named “Modified Neighbor\nGWAS for LASSO”, SNPs on strong LD (i.e., r2 > 0.8) were pruned before LASSO\nregression to reduce multicollinearity in the multiple regression by LASSO. We\ntherefore report the SNP position and MAF without LD for LASSO-selected SNPs. P-\nvalues are not applicable for LASSSO-selected SNPs as this is a method of model\nselection by machine learning.\nComment #15.\nLastly, it seems that table S4 may be missing. Could you please confirm this?\"\nResponse: We are sorry for our mistake and have now provided Supplementary\nData 1 that include a list of candidate genes in a .xlsx file format. In case this\nsupplementary file cannot be converted well in the submission system, the same\nresults are available as Table S4 via the bioRxiv preprint version at\nhttps://www.biorxiv.org/content/10.1101/2023.05.19.541411v1.supplementary-\nmaterial. We have also conducted a permutation test in response to comment #11\nbut none of the SNPs reported in Supplementary Data 1 became more significant\nthan the empirical significance threshold. We therefore avoid drawing strong\nconclusion from these candidate genes and have designated the neighbor genotype\neffects as polygenic as mentioned in the previous version.\nTo Reviewer #3\nComment #16.\nThe ms by Sato and colleagues is an extension of previously published work (Sato et\nal. 2021), where now data from an additional field site are analyzed. The ms is well\nwritten and the authors nicely demonstrate that effects of neighboring genotypes can\ninfluence traits in focal genotypes.\nThis is an important finding as it moves the quantitative analyses of genotype-\nphenotype relationships to a community level. The data of the two field sites also\nhighlight that these effects differ in differ environments, which is expected, but still\ngood to see.\nResponse: We are grateful to reviewer 3 for precisely evaluating the importance of\nour study. As evaluated here, our work provides “an important finding as it moves the\nquantitative analyses of genotype-phenotype relationships to a community level”.\nThe revised longer version now makes such a quantitative genetic and community-\nlevel significance clearer by accommodating the background and discussion in more\ndetail. Introduction and Discussion now include the following three revised\nparagraphs “Despite long-standing interest4,27, few studies have employed genome-\nwide polymorphism data to investigate stand-level properties in biodiversity research.\nwide association studies (GWASs) to dissect the genetic basis underlying stand-level\ngrowth within the model plant Arabidopsis thaliana28–30. For example, studies on 98\nA. thaliana genotypes have reported quantitative trait loci that mediate neighbor\ngenotype effects on plant growth30 and consequently mitigate belowground\ncompetition31. However, these studies require considerable effort in pairwise\ncultivation to sufficiently control neighbor composition among many genotypes. Due\nto this combinatorial cultivation effort, previous studies focused on limited pairs\nbetween ten focal and 98 counterpart genotypes in a controlled environment29,31.\nThese practical limitations remain an obstacle to enable large-scale GWAS and\nthereby identify the most beneficial pairs in anti-herbivore resistance in field\nenvironments.” (the third paragraph of Introduction, Line 65-76); “In this study, we\nsuccessfully predicted key genotype pairs that underpin associational resistance\nfrom numerous combinations of genotypes in which associational susceptibility was\nby far more prevalent. Intraspecific associational susceptibility to specialist\nherbivores is ubiquitous across plant species, as suggested by a previous meta-\nanalysis that reported the prevalence of negative effects of plant genotypic diversity\non resistance to specialist herbivores9. The present work overcame the difficulty in\ndistinguishing associational resistance from associational susceptibility by leveraging\ngenome-wide polymorphism data. Other than anti-herbivore resistance, recent\ngenome-wide studies have identified key genotype pairs that increase stand-level\ngrowth based on the limited number of genotype pairs in a controlled\nenvironment29,31. In contrast to this previous work, our approach can be applied to\nrandomized mixtures of many genotypes over a large space in the field. Such\napplicability to large-scale field experiments provides a novel solution to ecologically\nrelevant genetics and the prediction of intraspecific biodiversity effects on various\ntraits of interest.” (the first paragraph of Discussion, Line 275-288); and “Genomic\nprediction enabled us to identify key genotype pairs even when the target phenotype\nhad a polygenic basis. This strategy is comparable to genomic selection in plant\nbreeding, in which elite genotypes can be selected without identifying genes\nresponsible for a heritable but polygenic phenotype48–50. Identification of elite\ngenotypes is successful when genome-wide polymorphisms are sufficiently dense to\nrepresent the genetic potential of each genotype48. Prediction of genetic potential is\noften achieved using LASSO and its varieties49. Combining LASSO and Neighbor\nGWAS, our study expanded the concept of genomic selection towards neighbor\ngenotype effects and thereby showed the power of genome-wide polymorphism data\nto predict elite genotype pairs.” (the second paragraph of Discussion, Line 289-297).\nComment #17.\nI only have one major concern with the proposed model that I feel would be\nworthwhile to adders. One of the drawbacks of the analyses is that the\nneighborGWAS does not identify significant hits per se (see Fig. 2) and the authors\nargue that nonetheless they identify meaningful associations by combining the\nresults of the neighborGWAS with a selection scan. Here, a significant overlap of the\ntop SNPs identified with both methods is reported. As nice at this is, I feel it is rather\nindicative and I would prefer more direct evidence of associations.\nResponse: We agree that major findings should not be based on non-significant\nhits. Apart from the selection scan, we examined the number and frequency of\nestimated 𝛽* coefficients (histograms and boxplots in Figure 3) in more detail and\n\"\nfound that SNPs with negative 𝛽* (i.e., associational susceptibility) had higher MAF\n\"\nthan those with positive 𝛽* (associational resistance). As this pattern seems\n\"\nconsistent with the prevalence of genotype pairs with associational susceptibility at\npair level (Fig. 4a), we focused more on the result of GWAS itself (Results, Line 187-\n205). Instead, we moved the results of the selection scan from the main text to the\nsupplementary information (Supplementary Note 2 and Supplementary Figure 9).\nThe text in the revised Results section now reads “In the Neighbor GWAS, the sign\nof the estimated neighbor genotype effects 𝛽* represents positive or negative\n\"\ninteractions between the two alleles of paired neighbors, which corresponds to\nassociational resistance or susceptibility to herbivore damage and insect community\ncomposition34,43 (right inset in Fig. 1a; Supplementary Fig. 7c-f). We therefore\nexamined the number and frequency of SNPs with positive or negative 𝛽* separately\n\"\nfor the four phenotypes per site. We detected both negative and positive 𝛽* in the top\n\"\n0.1%-associated SNPs (histograms in Fig. 3e-h and m-p; Supplementary Fig. 8a-h).\nThe SNPs with associational resistance (𝛽* > 0) involved more minor alleles than\n\"\nthose with associational susceptibility (𝛽* < 0) in five out of the four phenotypes ×\n\"\ntwo sites (boxplots in Fig. 3e-h and m-p; Supplementary Fig. 8i-p), including the\nherbivore damage in Zurich (Fig. 3e). These data suggest that associational\nresistance may be less frequent among the A. thaliana genotypes. This frequency\nbias further motivated us to analyze another possible difference between the SNPs\nwith positive and negative 𝛽* , such as patterns of selection scanned by using\n\"\ngenome-wide polymorphism data. This selection scan revealed that SNPs with\nassociational resistance (𝛽* > 0) had more signatures of balancing selection than\n\"\nthose with associational susceptibility (𝛽* < 0) (Supplementary Fig. 9;\n\"\nSupplementary Note 2). These results highlight the different features of the top-\nscoring SNPs with associational resistance and susceptibility. Among many possible\nSNPs, we then attempted to determine key predictors that explain the herbivore\ndamage or insect community composition.” (Results, Line 187-205). As direct\nevidence of associations was difficult to obtain, we confirmed the significance of PVE\nusing permutation tests. Please see also our reply to the next comment.\nComment #18.\nOne problem to me is that it is not clear what a good null hypothesis for the\nneighborGWAS is. Hence significance thresholds are hard to estimate. A possible\nsolution could be permutations, where the neighborGWAS would be performed with\nthe same focal plant, but the neighboring (j=4 or 12) plants would just be randomly\nchosen instead of using the actual neighboring plants. These analyses should\nenable the calculation of a permutation based significance threshold instead of an\nsomewhat arbitrary enrichment store, as one would not expect an enrichment for the\npermutations.\nResponse: We appreciate this suggestion and have performed permutation tests in\nwhich neighboring plants are shuffled within field blocks. We conducted 99\npermutations for a subset of the present data and 20 permutations for the full data on\nherbivore damage in Zurich, because a single iteration takes approximately 13 hours\nper iteration per phenotype and thus is difficult to permute many times. We randomly\nshuffled neighboring plants 99 times as suggested and determined the 5% empirical\nthreshold of –log10(p) as 6.738. This is near to the Bonferroni threshold at p = 0.3,\nnamely –log10(p) = 6.783 in the case of 1,819,577 SNPs. While this empirical\nthreshold might not be equal among all the phenotypes, we could not find any\nsignificant SNPs above the threshold of –log10(p) at 6.78. When permuting\nneighboring plants, we also calculated PVE for the herbivore damage in Zurich and\nfound that observed PVE was above the 5% threshold of permuted ones. These\nresults of permutation tests confirmed the significant PVE without non-significant\nSNPs, supporting our conclusion about the polygenic basis of neighbor effects.\nWe have now included this additional result into Supplementary Note 1 and\nSupplementary Figure 6 as follows: “Supplementary Note 1. Permutation test for\nthe genome-wide significance threshold\nBonferroni correction is often considered conservative in GWAS because\nlinkage disequilibrium (LD) makes individual SNPs not statistically independent85,86.\nPermutation tests solve this issue by empirically determining the significant\nthresholds, although they are computationally demanding86. Using permutation tests,\nwe attempted to determine the significance level for the proportion of phenotypic\nvariation explained (PVE) by neighbor effects and GWAS. To test whether the\nobserved PVE was more significant than random, we randomly permuted\nneighboring plants within blocks and calculated the PVE for the herbivore damage in\nZurich at 𝐽 = 4. The net contribution of neighbor genotype effects was determined as\nPVE by the Neighbor GWAS model – the genomic heritability ℎ\" (see Methods\n“Phenotypic variation explained by focal and neighbor genotype effects”). We\ncalculated PVEs for 99 permutations and tested whether the observed PVE was\nhigher than the 95 percentile of the permuted PVEs. In contrast to PVE, the full data\nfor GWAS permutations required a computational time of approximately 13 h per\niteration per phenotype in a 24-thread parallel calculation. To deal with this difficulty,\nwe conducted a large and small set of permutations using the full data and subset\ndata of the herbivore damage in Zurich at the number of neighboring individuals 𝐽 =\n4. For the full data, we permuted neighboring plants within blocks 20 times and\nobtained the maximum –log (𝑝) among the 20 iterations. For the subset data, we\n10\nfocused on the data collected from the first two blocks (= 200 individuals × 2 blocks)\nin the Zurich site in 2017 because the highest SNP observed in these two blocks (–\nlog (𝑝) = 6.53) was very similar to that observed among all the phenotypes in this\n10\nstudy (–log (𝑝) = 6.55 for the internal feeders in Otsu at 𝐽 = 4; Fig. 3k;\nneighbor genotype effects 𝛽* was obtained per iteration. Their 95 percentile was\n\"\nconsidered an empirical significance level of 𝑝 = 0.05 for the subset permutation.\nAs a result of permutations, we found that the observed PVE by neighbor\ngenotype effects on the herbivore damage was higher than the 95 percentile of\npermuted PVEs in Zurich (Supplementary Fig. 6a), supporting the significance of\nPVE by neighbor genotype effects. When Neighbor GWAS at 𝐽 = 4 was applied to\npermuted full data, we found a higher association score of permuted –log (𝑝) than\n10\nthat of all tested SNPs for the herbivore damage in Zurich (max. –log (𝑝) = 7.60\n10\namong 20 permutations: Supplementary Fig. 6b). To further examine the genome-\nwide significance level of GWAS, we analyzed the subset data of herbivore damage\nobserved in Zurich with 99-times permutations of neighboring plants. This additional\npermutation on the subset data determined the 5% empirical threshold of –log (𝑝)\n10\nas ‘6.738’ at 𝐽 = 4 (Supplementary Fig. 6c). This empirical threshold is close to the\nBonferroni threshold at 𝑝 = 0.3, namely –log (𝑝) = ‘6.783’ in our case of 1,819,577\n10\nnone of the four phenotypes per site had more significant SNPs at 𝐽 = 4 than 6.783.\nThese supplementary results provide consistent evidence for significant PVE without\ngenotype effects.” Supplementary Figure 6 visualizes these results as follows.\nSupplementary Figure 6. Permutation tests for the genome-wide significance\nlevel and the proportion of phenotypic variation explained (PVE) by neighbor\ngenotype effects. a Histogram of 99-times permuted PVEs for the full data on\nherbivore damage in Zurich at 𝐽 = 4. PVE by neighbor genotype effects was\ndetermined as PVE by the neighbor GWAS model – the genomic heritability ℎ\",\nwhich showed the net contribution of neighbor genotype effects to phenotypic\nvariation. The vertical solid and dashed lines indicate the observed PVE and the 95\npercentile of the permuted PVEs, respectively. b 20-times permutation tests of the\nneighbor genotypic effects 𝛽 using the full data on the herbivore damage (leaf\n\"\nholes) in Zurich at 𝐽 = 4. The horizontal black line indicates the Bonferroni threshold\nat 𝑝 = 0.05. The horizontal gray line indicates the maximum –log (𝑝) among the\n10\npermutations. c 99-times permutation tests of the neighbor genotypic effects 𝛽\n\"\nusing the subset data on herbivore damage in Zurich at 𝐽 = 4. The horizontal black\nline indicates the Bonferroni threshold at 𝑝 = 0.05. The horizontal gray line indicates\nthe permuted threshold at 𝑝 = 0.05.\nEnd of the Letter\n19"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have addressed most of my comments but the point about classification of herbivores into external and internal feeders still does not make sense to me. As I pointed out in my previous review, normally herbivores such as leaf gallers, leaf miners and borers which feed inside plant tissues are called internal feeders, but the authors seem to use this term to refer largely to sucking insects (aphids and thrips) which feed externally. This is incorrect and this classification contradicts the herbivore classification used in plant-herbivore literature (see for instance recent review by Tooker and Giron https://www.frontiersin.org/journals/plant-science/articles/10.3389/fpls.2020.581816/full). Aphids and thrips are not internal feeders, they feed externally. Tooker and Giron discuss that most herbivorous insects in orders Thysanoptera and Hemiptera are external feeders (apart from few examples of gall inducers). Including them in the same group as leaf miners makes little ecological sense.\nIn addition, in sentence on line 463 the authors are talking about internal feeders (according to their definition) but call them external feeders: \"The individual number of external feeders to a large extent represented the individual number of sap-sucking herbivores (i.e., phloem-sucking aphids, cell content-sucking thrips, and sap sucking stink bugs)\". Division of herbivores into chewing and sucking would make more sense.\nMinor comment - line 321 - what is \"larval moths and butterflies\" - do you mean caterpillars?", "reviewer_2": "First of all, I appreciate the authors for their great effort in improving their manuscripts according to the comments of three reviewers. The manuscript seems to have improved and become clearer in terms of what they have done. On the other hand, I still don’t understand the novelty of this manuscript for publication in this journal. Although they claim that “Neighbor GWAS” is a new method (L78-L79), it was already published in 2021 by themselves. From their reply, the methodology has been established in their previous paper. This situation is very confusing. Instead of citing the previous study, the authors should clarify the novelty of their present study, even in the method section.\nL151-152: I suggested genome rotation instead of random sampling of genotypes from your accessions. This way, we can assess the effect of population structure. Random sampling typically introduces significance to the real test [1]. See “Genomic patterns of polymorphism” in the method section.\nL185-186: The neighbor effect might be polygenic, but there are many reasons why GWAS did not detect any significant peaks (allelic heterogeneity or false negative errors due to correlation with population structure, etc.,). How can the authors conclude it is a polygenic effect?\nL146-147: If I understand correctly, this is SNP-based heritability (from L563). Have the authors tried model selection with and without sigma_2^2? Neighbor genotype provides more genetic variation to explain the phenotypic variation uncaptured in the first sigma_2. Thus, an increased proportion of variance explained seems statistically expected when using the right genotypes.\nFig 2: The phenotypic distribution does not follow a normal distribution although all the models assume it (L504). Have you transformed the phenotypes to achieve the most normal distribution? (perhaps Box-cox transformation could be useful)\nFig 4b: I’m not sure if the correlation can be analyzed as a linear correlation. What is the reason that the authors conclude that the four data points on the right side are outliers? Also, the value of r = -0.37 seems inconsistent with the value in text L380.\nFig 4d: The authors selected several pairs to demonstrate the performance of their genetic prediction. Maybe I missed the information, but how many pairs have the authors tested? For example, has Bg-2 not shown a strong interaction with other lines as the authors predicted? If the authors have performed the test for other pairs, it should be included in the main figure to avoid cherry-picking. If not, only three pairs would not be sufficient to reach concrete conclusions.\n1. Nordborg M, Hu TT, Ishino Y, Jhaveri J, Toomajian C, Zheng H, et al. The pattern of polymorphism in Arabidopsis thaliana. PLoS Biol. 2005;3: e196.", "reviewer_3": "I appreciate the modifications made by the authors. To me, the ms greatly improved and my concerns were appropriately addressed. I don't have any further comments/suggestions."}, "author_response": "Italic: comments on the revised version; Grey italic: comments on the initial version;\nPlain: our reply; Blue: main text; blue underline: second revision; strike through: deleted\n----------------------------------------------------------------------------------------------------------------------\nReviewer #1 (Remarks to the Author):\nComment #1.\nThe authors have addressed most of my comments but the point about classification of\nherbivores into external and internal feeders still does not make sense to me. As I pointed\nout in my previous review, normally herbivores such as leaf gallers, leaf miners and borers\nwhich feed inside plant tissues are called internal feeders, but the authors seem to use this\nterm to refer largely to sucking insects (aphids and thrips) which feed externally. This is\nincorrect and this classification contradicts the herbivore classification used in plant-herbivore\nliterature (see for instance recent review by Tooker and Giron\nhttps://www.frontiersin.org/journals/plant-science/articles/10.3389/fpls.2020.581816/full).\nAphids and thrips are not internal feeders, they feed externally. Tooker and Giron discuss\nthat most herbivorous insects in orders Thysanoptera and Hemiptera are external feeders\n(apart from few examples of gall inducers). Including them in the same group as leaf miners\nmakes little ecological sense.\nResponse: Reviewer 1’s clarification is much appreciated. Now we see the point of\n‘exophages and endophages’ by reading Tooker and Giron (2020). We agree that our\nprevious wording is contradictory to the entomological terminology of external (i.e.,\nexophagous) and internal (endophagous) herbivores. To avoid entomologically incorrect\nterms, we have reworded the external and internal feeders in the previous version as\n“external chewers” and “other herbivores”, throughout the manuscript (Lines 116, 121, 164,\n166, 170, and 177; Main figures 2 and 3, Supplementary Figures, and Supplementary\nTables).\nOur primary aim was to classify herbivores into those involved in jasmonic acid (JA) and\nsalicylic acid (SA) pathways in plant defense responses, as suggested by Reviewer 1’s\nprevious comment #5 saying “Classifying herbivores into chewing and sucking would make\nmore sense as these herbivores are also known to trigger different plant defensive pathways\n(JA vs SA)”. Based on this JA-SA criterion, we attempted to classify all observed herbivores\ninto chewers and suckers, but leaf miners are known to elicit both JA and SA pathways (Yang\net al. 2021). It was also incorrect to call leaf minders non-chewers because they possess a\nmouthpart to chew plant tissues inside a leaf (reviewed by Tooker and Giron 2020). To keep\nour primary aim without using entomologically incorrect terms, we would therefore like to\nreword the internal and external feeders as ‘external chewers’ and ‘other herbivores’ in the\nrevised manuscript. Specifically, the revised Results and Methods section now reads “The\nfact that the major herbivore species observed in our study were specialists of Brassicaceae\n(Supplementary Fig. 1; Supplementary Table 2) led us to classify these insects into external\nchewers (i.e., herbivores that harbor outside and chew plant tissues) and other herbivores.\nThe former is known to trigger the defense pathways mediated primarily by the plant\nhormone jasmonic acid, and the latter include species that induce defense through salicylic\nacid (Zuest et al. 2016; Mertens et al. 2021; Yang et al. 2021).” (Results, Line 113-119); and\n1\n“The individual number of external chewers or other herbivores was defined as the total\nnumber of individuals of leaf-chewing species (e.g., beetles and caterpillars) or species that\nate internal parts of a plant (e.g., phloem-sucking aphids, cell-content-sucking thrips, sap-\nsucking stink bugs, and leaf miners). The reason for this classification was to reflect the\ndifference in plant defense responses to external chewers and other herbivores through\njasmonic acid and salicylic acid pathways, respectively (Zuest et al. 2016; Mertens et al.\n2021; Yang et al. 2021). Leaf miners are known as an endophagous herbivore (Tooker &\nGiron 2020) that can elicit both the jasmonic and salicylic acid pathways in plant defense\nresponses (Yang et al. 2021). Based on their relevance to plant defense responses mediated\nby the salicylic acid pathway, leaf miners were included in the category of other herbivores.\nThe individual number of other herbivores to a large extent represented the individual\nnumber of sap-sucking herbivores (i.e., phloem-sucking aphids, cell content-sucking thrips,\nand sap sucking stink bugs) because leaf miners were very rare (0.35% of all individuals in\nthe other herbivore category; see the raw data64).” (Methods, Line 463-476).\nRefs:\nTooker, J. F., and Giron, D. (2020) The evolution of endophagy in herbivorous insects.”\nFrontiers in Plant Science 11: 581816.\nYang, J. N., Wei, J. N., & Kang, L. (2021). Feeding of pea leafminer larvae simultaneously\nactivates jasmonic and salicylic acid pathways in plants to release a terpenoid for indirect\ndefense. Insect Science, 28(3), 811-824.\nComment #2.\nIn addition, in sentence on line 463 the authors are talking about internal feeders (according\nto their definition) but call them external feeders: “The individual number of external feeders\nto a large extent represented the individual number of sap-sucking herbivores (i.e., phloem-\nsucking aphids, cell content-sucking thrips, and sap sucking stink bugs)”. Division of\nherbivores into chewing and sucking would make more sense.\nResponse: We thank Reviewer #1 for pointing this out. This previous Line 463 included a\ntypo and should have been described as ‘The individual number of ‘internal’ feeders to a\nlarge extent represented…’. We have reworded ‘internal feeders’ as ‘other herbivores’, and\ncorrected this sentence as follows “The individual number of other herbivores to a large\nextent represented the individual number of sap-sucking herbivores (i.e., phloem-sucking\naphids, cell content-sucking thrips, and sap sucking stink bugs) because leaf miners were\nvery rare (0.35% of all individuals in the other herbivore category; see the raw data64).”\n(Methods, Line 472-476).\nComment #3.\nMinor comment - line 321 - what is \"larval moths and butterflies\" - do you mean caterpillars?\nResponse: Yes, we have reworded it as ‘caterpillars’ (Line 327).\n2\n----------------------------------------------------------------------------------------------------------------------\nReviewer #2 (Remarks to the Author):\nComment #4.\nFirst of all, I appreciate the authors for their great effort in improving their manuscripts\naccording to the comments of three reviewers. The manuscript seems to have improved and\nbecome clearer in terms of what they have done. On the other hand, I still don’t understand\nthe novelty of this manuscript for publication in this journal. Although they claim that\n“Neighbor GWAS” is a new method (L78-L79), it was already published in 2021 by\nthemselves. From their reply, the methodology has been established in their previous paper.\nThis situation is very confusing. Instead of citing the previous study, the authors should clarify\nthe novelty of their present study, even in the method section.\nResponse: We appreciate Reviewer #2’s further comments on our methodology. We agree\nthat ‘new method’ is misleading and have reworded it as ‘recently developed method’ (Line\n79 and 100). The present study is not a method paper but an application of the recently\ndeveloped method for empirical data. The novelty of our finding was described in the abstract\nas “Here, we identify key genotype pairs responsible for associational resistance to herbivory\nusing the genome-wide polymorphism data of the plant species Arabidopsis thaliana.” (Line\n32-34). Such novelty was also pointed out by the other reviewers. Reviewer #1’s initial\ncomment as follows “…, this study provides an interesting and promising proof of concept\nenabling to predict impacts of intraspecific mixed planting on herbivory and there is a\npotential to use such strategy in the future for crops.” Reviewer #3 recognized that the\npresent study was an extension of the previous work (Sato et al. 2021 Heredity) and\nevaluated its significance as follows “The ms by Sato and colleagues is an extension of\npreviously published work (Sato et al. 2021), where now data from an additional field site are\nanalyzed. The ms is well written and the authors nicely demonstrate that effects of\nneighboring genotypes can influence traits in focal genotypes.\nThis is an important finding as it moves the quantitative analyses of genotype-phenotype\nrelationships to a community level.” Therefore, the present study brings a new empirical finding\nbeyond the methodology alone.\nComment #5.\nL151-152: I suggested genome rotation instead of random sampling of genotypes from your\naccessions. This way, we can assess the effect of population structure. Random sampling\ntypically introduces significance to the real test [1]. See “Genomic patterns of polymorphism”\nin the method section.\n1. Nordborg M, Hu TT, Ishino Y, Jhaveri J, Toomajian C, Zheng H, et al. The pattern of\npolymorphism in Arabidopsis thaliana. PLoS Biol. 2005;3: e196.\nResponse: According to Reviewer #2’s initial comments, he/she did not mention ‘genome\nrotation’ but ‘permutation.’ Please see again Reviewer #2’s initial comment #11, which says\n“Additionally, have the authors conducted permutation tests for Neighbor GWAS? I believe\nit's essential to empirically validate \\beta_2 by random sampling to establish its distribution.”\ngenotypes in his/her initial comment #18, as follows “One problem to me is that it is not clear\nwhat a good null hypothesis for the neighborGWAS is. Hence significance thresholds are\nhard to estimate. A possible solution could be permutations, where the neighborGWAS would\nbe performed with the same focal plant, but the neighboring (j=4 or 12) plants would just be\nrandomly chosen instead of using the actual neighboring plants. These analyses should\nenable the calculation of a permutation based significance threshold instead of an somewhat\narbitrary enrichment store, as one would not expect an enrichment for the permutations.”\nTherefore, we conducted a permutation test following Reviewer #3’s concrete suggestion in\nthe first revision.\nIn response to this second comment #5 by Reviewer #2, we have conducted a permutation\ntest based on the genome rotation scheme (Supplementary Fig. 6d-f), but this did not change\nour conclusion. The genome rotation determined the genome-wide threshold above all SNPs\nfor the full set and subset data on the herbivore damage at Zurich, indicating no significant\nSNPs (p > 0.1: Supplementary Fig. 6e and f). The genome rotation-based PVE determined\nits 95% threshold at 0.045 (Supplementary Fig. 6d), confirming the significance of the\nobserved PVE of the herbivore damage. We have added these results to the revised Results,\nMethods and Supplementary Notes 1 as follows “Permutation tests confirmed that the\nobserved fraction of herbivore damage explained by neighbor genotype effects was greater\nthan that of randomly permuted neighboring plants in Zurich (permutation test, 𝑝 < 0.05;\nSupplementary Fig. 6a; Supplementary Note 1). This was also supported by a permutation\nscheme called genome rotation (Nordborg et al. 2005; Sasaki et al. 2022) that randomized\npopulation genetic structure (permutation test, 𝑝 < 0.05; Supplementary Fig. 6d;\nSupplementary Note 1).” (Results, Line 152-157); “This GWAS did not detect any significant\nSNPs for any of the four phenotypes at each site (𝑝 > 0.1 after Bonferroni correction: Fig. 3).\nThis was also supported by permutation tests that shuffled neighboring plants and by the\ngenome rotation scheme (𝑝 > 0.1; Supplementary Fig. 6b-c and e-f; Supplementary Note 1).”\n(Results, Line 186-189); “To empirically determine the significance threshold, we also\nperformed a permutation test in which neighboring plants within blocks or genomic positions\nwere randomized (Supplementary Fig. 6; Supplementary Note 1)” (Methods, Line 602 and\n623); “In addition to randomization of neighboring plants, we also employed a genome\nrotation scheme44,45 to randomize the population structure. This permutation scheme was\nconducted for PVEs and 𝑝-values using the subset and full data of the herbivore damage in\nthe same manner as described above.” (Supplementary Note 1); and “In addition to the\nrandomization of neighboring plants, we conducted permutation tests following a genome\nrotation scheme, which randomizes the population structure. This genome rotation-based\npermutation test also showed that the observed PVE for the herbivore damage was larger\nthan genome-rotated PVEs (red histogram and vertical line in Supplementary Fig. 6d) and\ndetermined the genome-wide significance threshold above all SNPs (red horizontal lines in\nSupplementary Fig. 6e and f).” (Supplementary Note 1). Supplementary Figure 6 shows\nthese supplementary results of the genome rotation (red lines in d-f) as follows:\nproportion of phenotypic variation explained (PVE) by neighbor genotype effects. a\nand d Histograms of 99-times permuted PVEs for the full data on herbivore damage in Zurich\nat 𝐽 = 4. PVE by neighbor genotype effects was determined as PVE by the neighbor GWAS\nmodel – the genomic heritability ℎ!, which showed the net contribution of neighbor genotype\neffects to phenotypic variation. The vertical solid line indicates observed PVE. a The blue\nhistogram and vertical dashed line represent the neighbor-permuted PVEs and their 95\npercentile. d The red histogram and vertical dashed line represent the genome-rotated PVEs\nand their 95 percentile. b and e 20-times permutation tests of the neighbor genotypic effects\n𝛽 using the full data on the herbivore damage (leaf holes) in Zurich at 𝐽 = 4. The horizontal\n!\nblack line indicates the Bonferroni threshold at 𝑝 = 0.05. The horizontal red and blue lines in\nb and e indicate the maximum –log (𝑝) of the neighbor- and genome-rotated 𝑝-values,\n10\nrespectively. c and f 99-times permutation tests of the neighbor genotypic effects 𝛽 using\n!\nthe subset data on herbivore damage in Zurich at 𝐽 = 4. The horizontal black line indicates\nthe Bonferroni threshold at 𝑝 = 0.05. The horizontal red and blue lines in c and f indicate the\n95 percentiles of the maximum –log (𝑝) among neighbor- and genome-rotated 𝑝-values,\n10\nrespectively.\nIn general, permutation tests are not common in GWAS because they impose a high\ncomputation cost. Because of this computational cost, many iterations are impossible unless\nspecialized hardware (e.g., high-performance multicore CPU or GPU) is available (John et al.\n2022). As far as we searched for other publications that employed the genome rotation as\nsuggested above (Nordborg et al. 2005), genome rotation (also called genome shuffling) was\nused for pre- or post-GWAS analyses, such as PVE estimates (Sasaki et al. 2015) and\nenrichment tests (Sasaki et al. 2022), but not for GWAS itself. We therefore think that our two\npermutation schemes using the full dataset and subset data are now sufficient to address the\nissue of empirical thresholds.\nRefs:\n(2022). Efficient permutation-based genome-wide association studies for normal and skewed\nphenotypic distributions. Bioinformatics, 38(Supplement_2), ii5-ii12.\nNordborg, M., Hu, T. T., Ishino Y., Jhaveri, J., Toomajian, C., Zheng, H., ... & Bergelson, J.\n(2005). The pattern of polymorphism in Arabidopsis thaliana. PLoS Biology, 3(7), e196.\nSasaki, E., Zhang, P., Atwell, S., Meng, D., & Nordborg, M. (2015). “Missing” G x E variation\ncontrols flowering time in Arabidopsis thaliana. PLoS Genetics, 11(10), e1005597.\nSasaki, E., Gunis, J., Reichardt-Gomez, I., Nizhynska, V., & Nordborg, M. (2022). Conditional\nGWAS of non-CG transposon methylation in Arabidopsis thaliana reveals major\npolymorphisms in five genes. PLoS Genetics, 18(9), e1010345.\nComment #6.\nL185-186: The neighbor effect might be polygenic, but there are many reasons why GWAS\ndid not detect any significant peaks (allelic heterogeneity or false negative errors due to\ncorrelation with population structure, etc.,). How can the authors conclude it is a polygenic\neffect?\nResponse: We meant the polygenic architecture by the result that there was significant PVE\nwithout significant GWAS peaks. Other authors also attributed complex traits to the polygenic\narchitecture when they detected significant SNP heritability but did not find significant SNPs\nabove the genome-wide threshold (e.g., human intelligence by Davies et al. 2011; human\nsuicidal behavior by Otsuka et al. 2019; sorghum root microbiome by Deng et al. 2021). To\navoid using the term ‘polygenic,’ we have rewritten them simply as “the lack of significant\nSNPs associated with the neighbor effects” or “no significant SNPs were detected” in this\nsecond revision (Lines 213 and 295), such that the word ‘polygenic’ is now used only once in\nthe main text (Line 191).\nRefs:\nDavies, G., Tenesa, A., Payton, A., Yang, J., Harris, S. E., Liewald, D., ... & Deary, I. J.\n(2011). Genome-wide association studies establish that human intelligence is highly heritable\nand polygenic. Molecular Psychiatry, 16(10), 996-1005.\nDeng, S., Caddell, D. F., Xu, G., Dahlen, L., Washington, L., Yang, J., & Coleman-Derr, D.\n(2021). Genome wide association study reveals plant loci controlling heritability of the\nrhizosphere microbiome. The ISME journal, 15(11), 3181-3194.\nOtsuka, I., Akiyama, M., Shirakawa, O., Okazaki, S., Momozawa, Y., Kamatani, Y., ... &\nHishimoto, A. (2019). Genome-wide association studies identify polygenic effects for\ncompleted suicide in the Japanese population. Neuropsychopharmacology, 44(12), 2119-\n2124.\nComment #7.\nL146-147: If I understand correctly, this is SNP-based heritability (from L563). Have the\nauthors tried model selection with and without sigma_2^2? Neighbor genotype provides\nmore genetic variation to explain the phenotypic variation uncaptured in the first sigma_2.\nThus, an increased proportion of variance explained seems statistically expected when using\nthe right genotypes.\nResponse: Yes, this is SNP heritability and we compared models with and without\nsigma_2^2 using the likelihood ratio test (see Line 578-583 in the version R1, and Line 463-\n465 in the initial submission). This model selection based on the likelihood ratio test aimed to\ntest whether the phenotypic variation uncaptured in the first sigma_2 was significantly\nexplained and PVE increased due to sigma_2^2, as Reviewer #2 pointed out here. To make\nthe point clearer, we have rewritten “genomic heritability” as “SNP heritability” (Line 577) and\ncorrected/revised the description of the likelihood ratio tests (Line 589-594) as follows:\nPrevious:\n1. Compute the null likelihood with 𝜎! = 0 and 𝜎! = 0.\n\" !\n2. Test the variance component 𝜎! in comparison with the null likelihood.\n\"\n3. Calculate ℎ! = 𝜎(!/(𝜎(!+𝜎(!).\n\" \" #\n4. Compute the focal likelihood with 𝜎! ≠ 0 and 𝜎! = 0.\n! !\n5. Test the variance component 𝜎! in comparison with the focal likelihood.\n!\n6. Calculate PVE = (𝜎(!+𝜎(!)/(𝜎(!+𝜎(!+𝜎(!).\n\" ! \" ! #\nRevised:\n1. Compute the likelihood with 𝜎! = 0 and 𝜎! = 0.\n\" !\n2. Compare models with and without 𝜎! based on the likelihood ratio test for 𝜎!.\n\" \"\n3. Calculate ℎ! = 𝜎(!/(𝜎(!+𝜎(!).\n\" \" #\n4. Compute the likelihood with 𝜎! ≠ 0 and 𝜎! = 0.\n\" !\n5. Compare models with and without 𝜎! based on the likelihood ratio test for 𝜎!.\n! !\n6. Calculate PVE = (𝜎(!+𝜎(!)/(𝜎(!+𝜎(!+𝜎(!).\n\" ! \" ! #\nComment #8.\nFig 2: The phenotypic distribution does not follow a normal distribution although all the\nmodels assume it (L504). Have you transformed the phenotypes to achieve the most normal\ndistribution? (perhaps Box-cox transformation could be useful)\nResponse: The phenotypic variation does not follow a normal distribution, whereas the\nmodel’s 𝑢 and 𝑒 assume normality. To clarify this point, we have revised the Line 504 as\nfollows.\nPrevious: In addition to the fixed effects, the random effects 𝑢 and residuals 𝑒 contribute to\n$ $\nphenotypic variation following a normal distribution as 𝑢 ∼ 𝑁(𝟎,𝜎!𝐊 +𝜎!𝐊 ) and 𝑒 ∼\n$ \" \" ! ! $\n𝑁(0,𝜎!).\n#\nRevised: In addition to the fixed effects, the random effects 𝑢 and residuals 𝑒 contribute to\n$ $\nphenotypic variation. These random effects and residuals follow a normal distribution as 𝑢 ∼\n$\n𝑁(𝟎,𝜎!𝐊 +𝜎!𝐊 ) and 𝑒 ∼ 𝑁(0,𝜎!). (Methods, Line 513).\n\" \" ! ! $ #\nquestion. There are two possible questions and answers.\nFirst, if his/her question is ‘whether’ we considered the transformation, our answer is yes. A\nnatural log-transformation was applied to all phenotypes before GWAS. This point was\ndescribed in both the initial and revised version (see Line 424-426 in the initial submission\nand Line 477-479 in version R1). Specifically, Line 487-489 in the Method section read “All\nphenotypes except for the leaf area loss (score variable) were ln(x+1)-transformed to\nimprove normality for GWAS and genomic prediction”. For further clarity, we have added the\nsame sentence to the caption of Figure 2.\nSecond, if his/her question is ‘why’ we used the log-transformation instead of Box-Cox\ntransformation, this was because the number of individual insects or holes was countable. In\nthe framework of generalized linear models (not generalized linear ‘mixed’ models), the\ncountable response should be analyzed using a Poisson or negative binomial error structure\n(e.g., Bolker et al. 2008; Faraway 2016), but these families are not available in the mixed-\nmodel framework in GWAS. Currently available mixed models for GWAS deal only with\nquantitative or binary traits (e.g., R packages named ‘milorGWAS’, ‘gaston’, and ‘COXMEG’).\nAs an alternative to the generalized mixed model framework, we used log-transformation as\na conventional option for count data analysis (e.g., Bartlett and Kendall 1946). To avoid zero\nwithin a logarithm, one was added to the countable phenotype as ln(x+1), because 1 is the\nminimum unit of count data and 0.5 individual does exist in the real world. In addition to\nGWAS, such a log transformation allows us to interpret LASSO-defined effect size (Line 701\nand Fig. 4a) as the log response ratio, i.e., log(a) – log (b) = log(a/b) between the mixture\nand monoculture.\nIn response to this second comment #8, we have checked no transformation and the Box-\nCox transformation and examined residuals of the linear mixed models. These residuals\nwere examined without fixed SNP effects because it was not feasible to check residuals for\nall GWAS SNPs (see “2.1.4 Transformation of quantitative phenotypes” in the\nSupplementary Information of Atwell et al. 2010 Nature). Please see the next page that shows\nresidual plots for the herbivore damage, the individual number of external chewers, the\nindividual number of the other herbivores, and the total number of insect species at Zurich\nand Otsu. These (non)-transformations made either normality or homoscedasticity worse.\nFor example, residuals of the individual number of herbivores worse deviated from normality\nthan those of log-transformed ones (Zurich l and r; Otsu l and r). Non-transformed herbivore\ndamage (leaf holes) at Zurich exhibited heteroscedasticity of residuals (Zurich a). The log-\ntransformation was therefore adopted because of the countable feature of response variables\nand no clear improvement by the other transformation methods.\n8\n9\n(reply to the comment #8 continued)\nMore generally, linear mixed models that assume a normal (Gaussian) error can be applied\nto any quantitative phenotype in GWAS. For example, 107 A. thaliana phenotypes analyzed\nby Atwell et al. (2010) included count data (e.g., the number of aphids). These phenotypes\nwere analyzed using a linear mixed model that assumed the Gaussian error, with no\ntransformation, log-transformation, or rank-transformation, which provided a reasonably\nuniform distribution of p-values (see “Transformation of quantitative phenotypes” in the\nSupplementary Information of Atwell et al. 2010). Therefore, residual normality should be\nimproved or made close to a normal distribution but this does not have to be completely\nnormal in GWAS (see Zhu et al. 2022 Methods in Molecular Biology; and Deng et al. 2022\nMethods in Molecular Biology for more recent GWAS protocols). Therefore, we think that our\nchoice of log-transformation is rational for the countable phenotypes.\nRefs:\nAtwell, S., Huang, Y. S., Vilhjálmsson, B. J., Willems, G., Horton, M., Li, Y., ... & Nordborg, M.\n(2010). Genome-wide association study of 107 phenotypes in Arabidopsis thaliana inbred\nlines. Nature, 465(7298), 627-631.\nBartlett, M. S., & Kendall, D. G. (1946). The statistical analysis of variance-heterogeneity and\nthe logarithmic transformation. Supplement to the Journal of the Royal Statistical Society,\n8(1), 128-138.\nBolker, B. M. (2008). Ecological models and data in R. Princeton University Press.\nDeng, S., Meier, M.A., Caddell, D., Yang, J., Coleman-Derr, D. (2022). Plant Microbiome-\nBased Genome-Wide Association Studies. In: Torkamaneh, D., Belzile, F. (eds) Genome-\nWide Association Studies. Methods in Molecular Biology, vol 2481. Humana, New York, NY.\nhttps://doi.org/10.1007/978-1-0716-2237-7_20\nFaraway, J. J. (2016). Extending the linear model with R: generalized linear, mixed effects\nand nonparametric regression models. Chapman and Hall/CRC.\nZhu, F., Fernie, A.R., Scossa, F. (2022). Preparation and Curation of Omics Data for\nGenome-Wide Association Studies. In: Torkamaneh, D., Belzile, F. (eds) Genome-Wide\nAssociation Studies. Methods in Molecular Biology, vol 2481. Humana, New York, NY.\nhttps://doi.org/10.1007/978-1-0716-2237-7_8\nComment #9.\nFig 4b: I’m not sure if the correlation can be analyzed as a linear correlation. What is the\nreason that the authors conclude that the four data points on the right side are outliers?\nAlso, the value of r = -0.37 seems inconsistent with the value in text L380.\nResponse: Both r values were correct. One value (r = –0.37) in the main text and Figure 4b\nshows the correlation coefficient when including four data points on the right side, whereas\nthe other value (r = –0.424) showed the coefficient when excluding the four data points. As\nwe were aware of the non-linearity, the four data points that appeared to deviate from\nlinearity were excluded to show the two cases of r values.\nfollowing sentence “Pearson’s correlation 𝑟 became larger (𝑟 = –0.424, *** 𝑝 < 0.001) when\nfour outlier genotypes were omitted from the x-axis at > 1” from Line 386 to avoid the\nconfusion. Instead of Pearson’s r that assume linearity and normality, we have also used a\nrank correlation, which yielded the almost same magnitude of correlation (rho = –0.38; see\nLine 251 and Fig. 4b).\nComment #10.\nFig 4d: The authors selected several pairs to demonstrate the performance of their genetic\nprediction. Maybe I missed the information, but how many pairs have the authors tested?\nFor example, has Bg-2 not shown a strong interaction with other lines as the authors predicted?\nIf the authors have performed the test for other pairs, it should be included in the main figure\nto avoid cherry-picking. If not, only three pairs would not be sufficient to reach concrete\nconclusions.\nResponse: Our major finding is the discovery of key genotype pairs that underpin\nassociational resistance to specialist herbivory (Results, Line 277-279). This finding itself is\nimportant because it was reported by a meta-analysis that intraspecific genetic diversity is\nless likely to mitigate specialist herbivory (Koricheva and Hayes 2018). We discussed this\npoint in the first revision as “Intraspecific associational susceptibility to specialist herbivores\nis ubiquitous across plant species, as suggested by a previous meta-analysis that reported\nthe prevalence of negative effects of plant genotypic diversity on resistance to specialist\nherbivores (Koricheva and Hayes 2018)” (Discussion, Line 283-286).\nWe did not plant any pairs other than the three pairs shown in the manuscript; therefore, our\nmixed planting experiment was not a cherry-picking that picked up three of many planted\npairs. To test the actual effects of mixed planting of the three pairs, it was necessary to\nmanipulate the mixture conditions and patch size simultaneously (Hambäck et al. 2014;\nUnderwood et al. 2014; Hahn and Orrock 2016). To do this, we planted 2,016 plants (1,008\nplants x 2 years; see Methods at Line 781), in which a single trial required an area of ca. 45\nm2. Therefore, cherry-picking itself was impossible in this large-scale field experiment, and\nwe think that we did our best in the fieldwork using the three pairs.\nRefs:\nHambäck, P. A., Inouye, B. D., Andersson, P., & Underwood, N. (2014). Effects of plant\nneighborhoods on plant–herbivore interactions: resource dilution and associational effects.\nEcology, 95(5), 1370-1383.\nHahn, P. G., & Orrock, J. L. (2016). Neighbor palatability generates associational effects by\naltering herbivore foraging behavior. Ecology, 97(8), 2103-2111.\nKoricheva, J, and Hayes, D. (2018) The relative Importance of plant intraspecific diversity in\nstructuring arthropod communities: a meta-analysis. Functional Ecology 32(7): 1704–17.\nUnderwood, N., Inouye, B. D., & Hambäck, P. A. (2014). A conceptual framework for\nassociational effects: when do neighbors matter and how would we know? The Quarterly\nReview of Biology, 89(1), 1-19.\n11\n----------------------------------------------------------------------------------------------------------------------\nReviewer #3 (Remarks to the Author):\nComment #11\nI appreciate the modifications made by the authors. To me, the ms greatly improved and my\nconcerns were appropriately addressed. I don't have any further comments/suggestions.\nResponse: We appreciate Reviewer #3’s helpful suggestions, which have greatly improved\nour manuscript.\n----------------------------------------------------------------------------------------------------------------------\nEnd of the Letter\nReviewer #2 (Remarks to the Author):\nThe authors have addressed all my concerns. Personally, I cannot provide positive comments on\nthis methodology (i.e., the genetic background noise is too high and not convinced two K matrices)\nand the significance of their discovery, but the editor and readers will evaluate the authors'\nfindings.\nAs my last comment on the authors' reply, 'In general, permutation tests are not common in\nGWAS because they impose a high computational cost,' maybe, it’s true, but the methodology\nperfectly fits to generate a null distribution of structural data. We cannot test our hypothesis\nwithout a clear null hypothesis statistically; thus, the test is incomplete without null disctirbution. I\nbelieve the authors know the permutation test is not suitable for testing structural populations.\nThe authors' effort in conducting genome-wide rotation sampling strengthened their conclusion.\nGood job!\nItalic: original comments; Plain: our reply\nREVIEWERS' COMMENTS\nReviewer #2 (Remarks to the Author):\nThe authors have addressed all my concerns. Personally, I cannot provide positive\ncomments on this methodology (i.e., the genetic background noise is too high and not\nconvinced two K matrices) and the significance of their discovery, but the editor and readers\nwill evaluate the authors' findings.\nAs my last comment on the authors' reply, 'In general, permutation tests are not common in\nGWAS because they impose a high computational cost,' maybe, it’s true, but the\nmethodology perfectly fits to generate a null distribution of structural data. We cannot test our\nhypothesis without a clear null hypothesis statistically; thus, the test is incomplete without null\ndisctirbution. I believe the authors know the permutation test is not suitable for testing\nstructural populations. The authors' effort in conducting genome-wide rotation sampling\nstrengthened their conclusion. Good job!\nResponse: We appreciate the last comment by Reviewer #2. Based on the three types of\ngenome-wide thresholds (i.e., the stringent Bonferroni threshold, the randomization of\nneighboring genotypes, and the genome rotation scheme), the present study has now\nreached the same conclusion that we did not detect any significant GWAS hits. We therefore\nthink that the risk of reporting false positive GWAS hits could be avoided in the present study.\nAs the Neighbor GWAS method is relatively new, we also recognize that its applicability\nshould be tested more in the future.\n1"}]

Decision letter