Disparate genetic variants associated with distinct components of cowpea resistance to the seed beetle Callosobruchus maculatus

Polygenic genome-wide association mapping identified two regions of the cowpea genome associated with different components of resistance to its major post-harvest pest, the seed beetle Callosobruchus maculatus. Cowpea (Vigna unguiculata) is an important grain and fodder crop in arid and semi-arid regions of Africa, Asia, and South America, where the cowpea seed beetle, Callosobruchus maculatus, is a serious post-harvest pest. Development of cultivars resistant to C. maculatus population growth in storage could increase grain yield and quality and reduce reliance on insecticides. Here, we use a MAGIC (multi-parent, advanced-generation intercross) population of cowpea consisting of 305 recombinant inbred lines (RILs) to identify genetic variants associated with resistance to seed beetles. Because inferences regarding the genetic basis of resistance may depend on the source of the pest or the assay protocol, we used two divergent geographic populations of C. maculatus and two complementary assays to measure several aspects of resistance. Using polygenic genome-wide association mapping models, we found that the cowpea RILs harbor substantial additive-genetic variation for most resistance measures. Variation in several components of resistance, including larval development time and survival, was largely explained by one or several linked loci on chromosome 5. A second region on chromosome 8 explained increased seed resistance via the induction of early-exiting larvae. Neither of these regions contained genes previously associated with resistance to insects that infest grain legumes. We found some evidence of gene–gene interactions affecting resistance, but epistasis did not contribute substantially to resistance variation in this mapping population. The combination of mostly high heritabilities and a relatively consistent and simple genetic architecture increases the feasibility of breeding for enhanced resistance to C. maculatus.


Introduction
Cowpea, Vigna unguiculata (L.) Walp., is a warm-season grain legume that serves as a major source of dietary protein, animal fodder, and soil fertility in arid and semi-arid parts of Africa, Asia, and South America (Ehlers and Hall 1997;Langyintuo et al. 2003;Asif et al. 2013). In these drought-prone regions, grain and fodder yields can be severely constrained by a variety of abiotic and biotic stresses, both in the field and in storage (Mishra et al. 2018;Boukar et al. 2020). The most important post-harvest pest of cowpea is the seed beetle Callosobruchus maculatus (F.) (Coleoptera: Chrysomelidae). Left unchecked, beetle populations can grow exponentially in storage, and cause up to 90% seed loss (Tarver et al. 2007;Deshpande et al. 2011). Infestations of C. maculatus are especially problematic for subsistence growers without adequate storage facilities or fumigants (Subramanyam and Hagstrum 2000;Mishra et al. 2019;Sanchez et al. 2019). A variety of non-insecticidal techniques have been developed and employed to control C. maculatus, including cultural, biological, and physical methods (Solleti et al. 1841;Cruz et al. 2016;Amusa et al. 2018). One approach has been to develop broadly resistant cultivars of cowpea and other grain legumes by first establishing the genomic basis of resistance (Schafleitner et al. 2016;Miesho et al. 2019). Candidate genes and QTLs conferring grainlegume resistance have often been identified by analyses of Communicated by Volker Hahn.
Identifying loci responsible for agronomically important traits can be enhanced by the use of so-called MAGIC populations, i.e., a series of crop lines developed from multi-parent, advanced-generation intercrosses (Cavanagh et al. 2008;Huang et al. 2015;Rollar et al. 2021). For cowpea, eight elite cultivars with desirable traits were subjected to structured crosses for eight generations, resulting in 305 homozygous recombinant inbred lines (RILs) (Muñoz-Amatriaín et al. 2017;Huynh et al. 2018, and references therein). Because each line carries a particular mosaic of genome blocks from the founding parents, RILs are highly suitable for distinguishing loci affecting agronomically relevant plant traits. The cowpea MAGIC population has already been used to identify QTLs associated with such traits as flowering time, plant growth rate, and seed size (Huynh et al. 2018). In this study, we use the F:8 cowpea RILs in a genome-wide association study to identify QTLs and candidate genes that may be involved in resistance to seed beetle damage.
A previous study determined that resistance to C. maculatus varied significantly among the eight parents used to establish the MAGIC population (Messina et al. 2019). Parental lines differed in their effects on larval survival and development time, egg-laying rates, and adult weight. Unexpectedly, a few cultivars also induced beetle larvae to exit seeds prematurely, which produced an added source of mortality because such larvae either fail to molt into viable adults or yield misshapen adults that do not mate or reproduce normally (Messina et al. 2019). Given this level of variation in seed-beetle resistance among the parents, we used similar assays here to measure variation in beetle performance in seeds of 293 RIL progeny. Our objective was to identify genomic regions and candidate genes associated with decreased beetle performance, which could then be manipulated for crop improvement in breeding programs (Douglas 2018;Kpoviessi et al. 2019;Grazziotin et al. 2020).
Inferences about genetic mechanisms underlying crop resistance to insects may depend on the particular assay used to measure pest performance, as well as the source of the pest population (Appleby and Credland 2003). This may be especially important for a cosmopolitan pest like C. maculatus, which consists of geographic populations ("biotypes") that are known to differ considerably in a suite of behavioral, physiological, and life-history traits, as well as in the range of crop species that they can infest (Messina 1991;Messina et al. 2020). We therefore used two different assays that were intended to provide complementary estimates of cowpea resistance, as well two divergent pest populations.
One assay focused on seed quality per se by measuring larval performance under optimal conditions, i.e., with no larval competition within seeds. A second assay estimated pest performance and F1 progeny production in a more realistic, competitive environment. For the second assay, we included C. maculatus populations known to differ in traits affecting host use and rates of infestation Messina et al. 2020).

Development of the cowpea RILs
A detailed description of the formation of the MAGIC population can be found in Huynh et al. (2018). Briefly, eight parental cultivars (SuVita 2, CB27, IT93K-503-1, IT89KD-288, IT84S-2049, IT82E-18, IT00K-1263, and IT84S-2246 were chosen because they were high-yielding and possessed resistance or tolerance to various abiotic and biotic stresses (Huynh et al. 2018, and references therein). The eight parents also harbored considerable genetic diversity; genotyping indicated that about 68% of 51,128 putative single-nucleotide polymorphisms (SNPs) were polymorphic among the eight parents, and 11,848 SNPs were unique to a particular parent (Muchero et al. 2009(Muchero et al. , 2013Huynh et al. 2018). Parental cultivars were systematically inter-crossed for eight generations to yield 305 RILs (Muñoz-Amatriaín et al. 2017), each of which carried a mosaic of genome blocks from all founders. Because of transgressive segregation in their recombined genomes (de los Reyes 2019), the RILs exhibited wide phenotypic variation, in many cases beyond that observed among the eight parents (Huynh et al. 2018). To detect the possibility of similar transgression for resistance to seed beetles, we included in our experiments 293 RILs along with the eight parents, for a total of 301 cowpea genotypes.

Beetle life history and source populations
Eggs of C. maculatus are attached to the surfaces of dried seeds, and the hatching larva chews into the seed directly beneath the egg-laying site. All larval development must be completed within the single, natal seed. As a consequence, there is strong competition when multiple larvae must share a seed (Mitchell 1975;Messina 1991). Emerging adults commence mating and oviposition within hours after emergence. They require neither food or water and do not usually receive either in human stores of grain legumes. It is likely that C. maculatus was originally associated with the wild progenitor of cowpea in Africa (Huynh et al. 2013;Oas et al. 2015). Global transport of infested grain legumes (pulses) has introduced the beetle to all regions where grain legumes are stored, including areas where environments outside of storage facilities are unsuitable for maintaining beetle populations (e.g., Konrádhsdóttir et al. 2021). Because this pest has probably been associated with human stores of pulses for thousands of years (Tuda et al. 2006(Tuda et al. , 2014, it is especially well suited to assays in typical laboratory environments (Messina 2004).
We examined RIL suitability for seed beetles from two populations that were chosen because of their divergent life histories (Mitchell 1990;Messina 2004). The first two experiments used a beetle population derived from an infestation of mung bean (V. radiata [L.] Wilczek) and the related black gram (Vigna mungo [L.] Hepper) in southern India (hereafter = the SI population) (Mitchell 1991). It had been reared on mung bean (Berken cultivar) for > 350 generations in our laboratory at the start of the current experiment. Despite its long tenure in the laboratory, the SI population maintains genetic variation for a variety of host-use traits (Fox et al. 2009;Messina et al. 2009;Gompert and Messina 2016). Because of its chronic association with a smallseeded host (mung bean), SI beetles exhibit two unusual traits that could affect rates of population growth and seed infestation (Messina 2004). First, females are especially reluctant to place additional eggs on seeds already bearing a few eggs and thus produce a mostly uniform distribution of eggs among seeds (Messina and Mitchell 1989). If all available seeds bear two or three eggs, females will die without depositing all of their remaining eggs. Second, and more importantly, larvae exhibit especially strong contesttype competition within seeds, so that even larger host seeds receiving multiple eggs and larvae consistently yield only one or two adults (Messina 1991;Toquenaga and Fujii 1991;Toquenaga 1993;Fox and Messina 2018).
The third experiment instead used a population derived from infested cowpea in California, USA ( hereafter = the CA population ). It had been maintained on cowpea (California black-eyed pea cultivar) in the laboratory for > 150 generations at the start of the current study (Dowling et al. 2007;Tuda et al. 2014;Downey et al. 2015). Life-history and behavioral traits in the CA population are more typical for those of other geographic populations of C. maculatus. In contrast to SI females, CA females continue to add many eggs per seed (sometimes > 10 ) when there are few available seeds (Messina et al. 2020). Moreover, larvae exhibit a tolerant or scramble-type competition within seeds (e.g., Mano and Toquenaga 2008), so that a large, heavily infested seed can yield > 10 adults, albeit of relatively small body size.
We maintained beetle populations and conducted all experiments in a growth chamber at 25 • C and constant light. New generations of stock cultures were formed by adding 1500-2500 adults (estimated by volume) to a two-liter jar containing approximately 750 g of seeds, which corresponds to about 12,000 mung beans (for SI) or 3200 cowpeas (for CA). All seeds used to maintain stock cultures were organically grown and obtained in bulk quantities (about 55 kg per lot) from Azure Standard (Dufur, Oregon, USA).

Host suitability under optimal conditions
In the first experiment, we presented seeds to SI females in a randomized design and estimated larval performance on the 293 RILs and the eight parent lines. In each of 301 numbered Petri dishes, we added 10 randomly chosen seeds and a single pair of newly emerged beetles. Each seed within a dish was numbered with a permanent marker, and the combination of dish number and seed number indicated RIL identity. By adding a random assortment of seeds to each dish, we could reduce the effects of variation among beetle parents on estimates of larval performance and RIL suitability. Newly emerged beetle pairs were obtained by first sieving all adults in an actively emerging SI culture and then returning to collect beetles that had emerged within one hour after sieving. Thus, test females had not yet commenced oviposition when they were added to the dishes.
After 24 h, we inspected the dishes, removed all seeds that bore ≥ 1 egg, and placed those seeds into a second, labeled dish. The remaining seeds were exposed to females for an additional 24-h period. After this second oviposition period, we inspected all seeds in each of the two cohorts and used a metal probe to kill excess eggs so as to leave only one egg per seed. This step ensured that each larva developed with no competition within a seed. We discarded seeds that did not receive an egg during either exposure period (about 15% of the 3010 seeds). After 10-15 days, seeds bearing a single hatched egg, i.e., a single larva within the seed, were isolated in 4-ml, labeled glass vials (Messina et al. 2019).
Vials were inspected daily, and each newly emerged adult was sexed and weighed on an electronic balance (Mettler AE 160). We also recorded instances when the isolated seed yielded a prematurely exiting larva instead of a normal adult, as had been observed on three of the eight parental lines (Messina et al. 2019). Such larvae may die in the final larval stage or may molt outside the seed into a pupa or a misshapen, nonviable adult. We continued to inspect vials until 10 days after the last individual emerged. The three performance variables in this experiment were larval survival to adult emergence, development time from oviposition to adult emergence, and weight at adult emergence.

Beetle productivity under competitive conditions
Two additional experiments used the SI and CA populations to estimate larval performance and production of beetle adults in a competitive environment. As described above, these populations differ substantially in their responses to competition, in terms of both the tendency to add eggs to egg-laden seeds and the way larvae interact within seeds (contest vs. scramble competition). We therefore employed a protocol that would cause some degree of competition for both egg-laying sites and resources within seeds, as would typically occur as populations grow in untreated storage infestations.
In these experiments, each of the 301 dishes received 10 seeds of the same RIL or parent genotype instead of a random combination. We added two pairs of newly emerged beetles from either the SI (second experiment) or CA (third experiment) populations to each dish. Pairs were removed after three days. After a subsequent six days, when eggs had hatched, we counted the number of eggs laid by each pair of females per dish, as an indicator of host attractiveness for oviposition. Accurate egg counts can be obtained even after egg hatch because egg covers remain completely intact on the seed surface after larvae burrow into the seed immediately below the egg. A few dishes in each experiment received no eggs during the three-day exposure period: three dishes in the SI experiment and 13 dishes in the CA experiment (including a dish with a parental genotype, IT84S-2246). We excluded these dishes from the analyses of egg number because we could not distinguish whether the RILs were highly deterrent to egg-laying or received no eggs because of mating failure. The latter explanation is more likely because none of the 301 genotypes was so deterrent to egg-laying so as to receive no eggs in both experiments. In addition, two dishes in the CA experiment inadvertently did not receive beetles. Thus, we assayed 298 cowpea genotypes in the SI-population experiment and 286 genotypes in the CA-population experiment. Dishes were inspected daily, and we recorded the number of adult emergers or earlyexiting larvae.
Dependent variables included the number of eggs laid, survival to adult emergence (number of emerged adults/ number of eggs), and development time from oviposition to adult emergence. We standardized estimates of development time by assuming each egg was laid in the middle of the three-day oviposition period. The total number of emerged adults per dish (i.e., per RIL) was included as an additional, composite variable that depended on both the number of eggs laid and the proportion of larvae surviving to adulthood. Whereas survival in the first experiment (with no larval competition) depended solely on seed quality, survival in the second and third experiments could also depend on the severity of larval competition. For example, survival might be higher in a RIL that received relatively few eggs than in a RIL that was equally suitable for larvae but received many eggs. Nevertheless, by integrating host attractiveness and suitability for larvae, the number of emerged adults serves as an overall indicator of host susceptibility in a realistic scenario.

Polygenic genome-wide association mapping
We analyzed 32,130 polymorphic SNPs from the 293 RILs to identify genomic regions and candidate genes associated with the resistance traits from each experiment. Specifically, we fit Bayesian sparse linear mixed models (BSLMMs) with gemma (version 0.95alpha) to estimate the genetic contribution to each resistance trait (Zhou et al. 2013). Unlike traditional genome-wide association (GWA) mapping methods that test each genetic marker separately, the BSLMM method fits all SNPs in a single model and thus mostly avoids issues related to testing large numbers of null hypotheses. With this approach, trait values are determined by a polygenic term and a vector of the (possible) measurable effects of each SNP on the trait ( ) (Zhou et al. 2013). Bayesian Markov chain Monte Carlo (MCMC) with variable selection is used to infer the posterior inclusion probability (PIP) for each SNP, that is, the probability that each SNP has a nonzero effect, and the effect conditional on it being nonzero (Guan and Stephens 2011). The polygenic term denotes each individual's expected deviation from the mean phenotype based on all of the SNPs. This term accounts for phenotypic covariances among individuals caused by their relatedness or overall genetic similarity (Zhou et al. 2013). The kinship matrix also serves to control for population structure and relatedness when estimating effects of individual SNPs ( ) along with their PIPs. Similarly, SNPs in linkage disequilibrium (LD) with the same causal variant effectively account for each other, such that only one or the other is needed in the model, and this redundancy is captured by the posterior inclusion probabilities.
The hierarchical structure of the model makes it possible to estimate additional parameters that describe aspects of a trait's genetic architecture (Guan and Stephens 2011;Zhou et al. 2013;Lucas et al. 2018;Gompert et al. 2019). These include the percentage of the phenotypic variance explained (PVE) by additive genetic effects (which includes and the polygenic term, and should approach the narrow-sense heritability), the percentage of the PVE due to SNPs with measurable effects or associations (PGE, the percentage of the phenotypic variance explained by genic effects, which is based only on ), and the number of SNPs with measurable associations (n-). All of these metrics use MCMC to integrate over uncertainty in the effects of individual SNPs, including whether these are nonzero. Lastly, using this BSLMM approach, it is also possible to obtain genomicestimated breeding values (GEBVs), that is, the expected trait value for an individual from the additive effects of their genes, as captured by both and the polygenic term (Lucas et al. 2018;Gompert et al. 2019).
We used gemma to fit BSLMMs for the following resistance traits: (i) percent survival, development time, male weight and female weight for SI beetles under optimal conditions, (ii) percent survival, development time, total number of eggs, total number of adults, percent early-exiting larvae for SI in a competitive environment, and (iii) percent survival, development time, total number of eggs, total number of adults, percent early-exiting larvae for CA in a competitive environment. Traits were standardized prior to analysis, such that each had a mean of zero and standard deviation of one. For each resistance trait, we based our GWA parameter estimates on 10 MCMC runs, each comprising 1 million iterations and a 200,000 iteration burn-in. Every 10th MCMC sample was retained to form the posterior distribution. Genomic-estimated breeding values (i.e., polygenic scores) were then calculated from the parameter estimates and were used to calculate genetic covariances among traits.

Tests for epistasis
Our modeling approach (BSLMM) considers sets of genetic loci simultaneously, but would still fail to capture possible contributions of gene-gene interactions (i.e., epistasis) to resistance-trait variation. We therefore conducted complementary analyses to test for epistasis affecting these traits using the program MAPIT (Crawford et al. 2017). This method avoids the large combinatorial burden of testing for pairwise or higher-order epistasis; it instead evaluates the null hypothesis that the variance component for the vector of epistatic terms involving interactions between each SNP and all other SNPs is zero (Crawford et al. 2017). We computed p-values for tests of marginal epistasis using the recommended hybrid method that first implements a z-test to compute a p-value and then re-computes the p-value with the Davies method if the initial values is less than 0.05. (This step enhances the precision of calculations for low p-values without adding a large computational burden.) We then focused on the set of SNPs showing evidence of marginal epistasis: 139 SNPs with false-discovery rate q-values < 0.01 for development time (the trait with the most evidence of epistasis; see Results). We re-fit BSLMMs for development time while allowing for pairwise epistasis between all pairs of these SNPs (as in Nosil et al. 2020). This was done by including the products of the centered genotypes for each pair of SNPs as independent variables in the model. We were interested in whether including these terms increased the amount of trait variance explained by genetics, and whether the epistatic terms had non-trivial PIPs in the final models. We again based our GWA parameter estimates on 10 MCMC runs, each comprising 1 million iterations and a 200,000 iteration burn-in.

RIL suitability under optimal conditions
Larval performance varied considerably among the 301 cowpea genotypes. In the absence of competition, the average survival of SI larvae was high (80.5%), but it fell below 50% in about a tenth of the RILs and in two parental genotypes (Table 1, Fig. 1a). Most RILs conferred relatively rapid larval growth, with a mean egg-to-adult development time ≤ 25 days, but development required > 30 days for a substantial fraction of the RILs and for two parental genotypes (Table 1, Fig. 1b). This nearly disjunct distribution of mean development time was also observed in the pilot study of the eight parents. Similarly, for both males and females, there was an approximately

MAGIC RILs Parents
twofold difference in mean adult weight between the most and least suitable RILs (Fig. 1c,d). In most RILs, no larvae exhibited the aforementioned "suicidal" behavior by exiting seeds prematurely. However, 10-100% of larvae did so in two parental genotypes and about a fifth of the RILs (Table 1, Fig. 1e). In addition to the observed variation in RIL suitability for individual performance traits, phenotypic and genetic correlations indicated that many RILs could be classified as relatively good or poor hosts for multiple traits simultaneously, e.g., some conferred both poor survival to adult emergence and slow development. These correlations are elaborated below.

RIL suitability in a competitive environment
As expected, the divergent life-histories of the SI and CA populations produced stark differences in overall beetle performance under a competitive regime. When pairs of females were provided only 10 seeds for three days, CA females laid an average of 44.1 eggs, whereas SI females laid only 31.7 eggs. Despite lower subsequent densities of SI larvae compared to CA larvae, contest competition within seeds caused only 30.1% of SI larvae to survive to adult emergence. As a consequence, seeds yielded an average of only 9.1 emerging adults per 10 seeds, i.e., around one adult per seed, as expected from pure contest competition. In contrast, the tolerant behavior of CA larvae permitted an average survival to adult emergence of 78.6%, which was very similar to the average survival of SI larvae in the absence of competition. Moreover, the average number of emerging CA adults was nearly four times higher: 34.6 adults per 10 seeds. Unlike larval survival and the number of emerging adults, mean development time was similar between beetle populations (27.3 days for SI vs. 28.1 for CA).
The magnitude of variation in RIL suitability (or resistance) can be further illustrated by phenotypic correlations between populations (Fig. 2). For development time, a high positive correlation indicated that the two populations responded very similarly to variation among cowpea genotypes. At the other extreme, the number of eggs laid per RIL was not correlated between populations, i.e., RILs that elicited relatively greater egg-laying by SI females did not consistently do so among CA females, and vice versa. There were moderate positive correlations between populations for larval survival and for the number of emerging adults (Fig. 2a,d). Most RILs again yielded no early-exiting larvae in both the CA and SI populations, so that this trait was only weakly correlated between populations. As with the no-competition experiment, we describe below genetic correlations between individual traits, as well as between populations for a given trait (Fig. 2).

Overall genetic architecture
We detected segregating genetic variation affecting each of the resistance traits in the cowpea RILs (Table 2, Fig. 3). Bayesian point estimates of the percentage of trait variation attributable to additive genetic effects (PVE) ranged from 5% for the number of eggs laid to 64-67% for development time (both in the SI population). Under most conditions, development time and survival were highly heritable and were affected by a moderate number of genetic variants with measurable effects, i.e., the traits exhibited a high PGE with n = 2 -5 causal variants (Table 2, Fig. 3). In contrast, other traits, such as the total number of eggs or adults produced by the SI population, were less heritable and better explained by polygenic background effects. The heritability (PVE) of early-exiting ("suicidal") larvae depended on the C. maculatus population as well as experimental conditions; it ranged from 24-62%, but was consistently associated with a modest number of causal variants (3-5).
Estimates of genetic correlations among traits suggested similar genetic bases for some traits across populations and between competition regimes, but not for others. For example, the genetic correlation between survival of SI larvae with competition vs. without competition was 0.28. Genetic correlations between SI survival and CA survival were relatively high: 0.44 when SI larvae developed without competition, and 0.48 when SI larvae, like CA larvae, developed with competition (Fig. 4). However, the genetic correlation between the number of eggs laid and the subsequent number of F1 adults differed between populations; it was much higher in the "tolerant" CA population ( r = 0.85 ) than in the "contest" SI population ( r = 0.43).
We also detected moderately high genetic correlations between other metrics of cowpea resistance. For example, development time and survival consistently exhibited substantial negative genetic correlations, that is, RILs that conferred slow development also caused lower survival rates (for SI under no competition r = −0.36 , for SI with competition r = −0.57 , and for CA with competition r = −0.79 ). Similarly, development time was negatively correlated with the adult weight of both sexes when SI beetles developed without competition (the only experiment in which we measured weight at adult emergence) (for males r = −0.46 , for females r = −0.94 ). Thus, despite taking longer to emerge, slower-developing larvae reached a smaller adult size. In contrast, some components of resistance had largely independent genetic bases, as indicated by near-zero genetic correlations. Perhaps most notably, the proportion of earlyexiting larvae was mostly unrelated to survival or development time, especially in the two experiments that used a competitive regime: r = 0.02 for survival and −0.02 for 1 3 development time in the SI population, and −0.03 and −0.20 respectively in the CA population (Fig. 4).

Genotype-phenotype associations
We detected credible associations between individual SNPs and trait values for over half of the cowpea resistance traits (Figs. 5, 6). The most credible evidence for genotype-phenotype associations involved SNPs on cowpea chromosomes 5 and 8, with some weaker and more idiosyncratic evidence of associations on other chromosomes. (Here, we follow the chromosome numbering standard devised by Lonardi et al. 2019b.) For example, we observed credible associations (i.e., PIPs > 0.1 ) with SNPs on chromosome 5 for development time in all experiments, for survival in both the SI and CA populations in the presence of competition, for SI female weight in the absence of competition, and for the number of SI adults produced with competition (Figs. 5, 6).

MAGIC RILs Parents
In most cases, SNPs associated with these resistance traits mapped to approximately the same region of the genome, specifically between map positions 13.76 and 18.54 cM. This region includes 243 SNPs and spans 1.   5, the SNPs most associated with resistance traits (i.e., with PIPs > 0.2 ) resided within 10 distinct genes, including an ethylene-responsive element binding factor, aspartyl protease, and alpha/beta-hydrolase (Table 3). All but one of these genes are known to be expressed in cowpea seeds. For each of these 10 SNPs, we designated the allele associated with increased resistance based on the sign (i.e., direction) of the inferred phenotypic effect from the BSLMM. The full set of alleles conferring increased resistance at these loci was identified in haplotypes from three of the MAGIC RIL parents, IT89KD-288, IT84S-2246 and IT93K-503-1. In contrast, most SNPs associated with early-exiting larvae were on chromosome 8 (Fig. 5c,f,   one genetic variant causing the suicidal-larvae trait was between 0.47 (for the CA population) and 0.86 (for the SI population without competition). Notably, this genomic interval also contains several SNPs credibly associated with the total number of adults produced with competition in the CA population (Fig. 6). Within this region, the SNPs most associated with the early-exiting larvae resistance trait (PIPs > 0.2 ) resided within three genes, all of unknown function but known to be expressed in cowpea seeds at > 2 transcripts per million (Vigun08g171700, Vigun08g179000, Vigun08g180500). In this case, the full set of alleles conferring increased resistance at these loci was identified in haplotypes from two of the MAGIC RIL parents, IT84S-2246 and IT93K-503-1, both of which also contained resistance alleles on chromosome 5.

Evidence for epistasis
For most traits, we detected little compelling evidence of marginal epistasis (Fig. 7). Evidence of epistasis mostly fell short of genome-wide significance (i.e., p > 0.05 32130 ). We failed to obtain reliable results from tests of epistasis affecting early-exiting larvae; there was an excess of small p-values, likely stemming from the highly skewed distribution of these data. However, we did obtain significant evidence of marginal epistasis for six SNPs on chromosome 5 for development time in the CA population (Fig. 7f). These SNPs were at map positions 12.79 cM (one SNP) and 13.76 cM (five SNPs), and thus generally coincided with the region of the genome associated with additive-genetic variation for this trait (Fig. 5h). Moreover, a similar albeit weaker signal of marginal epistasis was observed in this same region for development time in the SI population (Fig. 7b,d).
We therefore focused further tests of epistasis on development time alone. Given the evidence of epistasis affecting this trait, we first converted p-values to falsediscovery rate q-values using the qvalue package (version 2.15.0) in R (version 3.6.3) (Storey et al. 2017). We then re-ran the polygenic BSLMM for each development time trait with interaction (epistasis) terms for all pairs of 139 SNPs with q-values less than 0.01 ( = 9591 interaction terms). The percent of trait variation explained by these new models with pairwise epistasis was similar to that of models without epistasis: 65% (90% ETPI 60-71%) with epistasis vs. 67% without epistasis for SI with no competition, 59% (90% ETPI 52-66%) vs. 64% without epistasis for SI with competition, and 58% (90% ETPI 50-64%) vs. 63% without epistasis for CA with competition (see also Table 2). Moreover, posterior inclusion probabilities for the epistasis terms never exceeded 0.001. Inclusion of epistasis terms thus did not improve model explanatory power even for development time.

Discussion
The objective of this study was to use a MAGIC cowpea population to identify the genomic basis of resistance to a highly destructive, storage pest (Huynh et al. 2018;Messina et al. 2019). In three large-scale experiments, we measured multiple aspects of seed-beetle performance on 286-301 F:8 genotypes. In so doing, we could examine the degree to which conclusions about crop resistance depended on experimental protocol, the particular pest trait measured, or the geographic source of the pest population. Few studies have considered each of these factors in attempts to determine the genomic basis of resistance to herbivorous insects. Our results revealed instances where genetic mechanisms of resistance appear to be general and robust, as well as cases in which conclusions did vary according to the particular experimental conditions. Bayesian models indicated that the amount of variation due to additive-genetic effects was relatively high for most but not all traits, and there was also trait-specific variation in the apparent location and number of causal genetic variants.

Genomic basis of cowpea resistance
Perhaps the clearest genetic signal was associated with eggto-adult development time and survival, two traits that were highly genetically correlated and would have a relatively large influence on the rate of beetle population growth in storage. These traits were also correlated with weight at adult emergence, which in turn is known to influence two additional fitness components: adult longevity and female fecundity (Messina 1991;Fox et al. 2004). Notably, SNPs associated with development time, survival, and weight all generally mapped to the same region on chromosome 5, with relatively few causal variants. Further research is needed to identify which of the genes that reside within this region of chromosome 5 are responsible for reducing the beetle performance (as discussed further below). It seems likely that these genes pleiotropically affect multiple beetle performance traits and hence the overall level of seed resistance. A second distinct genomic region conferring resistance occurred on chromosome 8, which contained most SNPs associated with the frequency of early-exiting larvae. This trait had not been recognized or included in previous studies of resistance to C. maculatus, but could represent an important, additional source of larval mortality (Messina et al. 2019). The independent genetic basis of this resistance trait was also confirmed by near-zero genetic correlations between the proportion of early exits on a RIL and either larval development time or survival. In a high-quality legume host, prepupal larvae ensure that there is only a thin layer of seed coat between their burrow and the outside of the seed. Newly molted adults only need to push open a piece of seed coat to emerge, and leave a smooth, circular exit hole (Southgate 1979). Early-exiting larvae instead burrow through the thin seed coat before pupation, and leave a distinct, jagged opening. It remains unclear which chemical or physical seed properties account for the aberrant behavior, but it was observed in all three experiments in this study, and was previously shown to cause larval mortality in the parents of the MAGIC population (Messina et al. 2019). Larvae in the earlier study developed without competition and the frequency of early exits was essentially bimodal; no larvae exited seeds in five parental lines, but a non-trivial fraction (20-36%) did so in three of the parents. Because inducing early exits was here shown to be heritable, the broadest and most durable level of cowpea resistance might be obtained by combining alleles simultaneously conferring slow development, low survival, small adult weight, and a non-trivial amount of pest mortality from early-exiting behavior.
Mapping analyses indicated that SNPs in ten annotated genes on chromosome 5 and three genes of unknown function on chromosome 8 were associated with cowpea resistance. The specific alleles conferring resistance were inherited from one or more of three closely related cultivars: IT89KD-288, IT84S-2246, and IT93K-503-1 (Huynh et al. 2018). We previously discovered that C. maculatus larval development time increased in each of these three cultivars, and development in IT93K-503-1 in particular was associated with lower survival and a high frequency of earlexiting larvae (Messina et al. 2019). None of these candidate genes had been identified in earlier studies of resistance to C. maculatus in cowpea or in other grain legumes. For example, Miesho et al. (2019) mapped cowpea resistance to C. maculatus in a panel of 217 cowpea accessions considered to be representative of worldwide cowpea diversity (Munoz-Amatriain et al. 2016). They identified six candidate genes associated with egg number, development time and insect emergence. None of the genes were on chromosome 5, where we found the strongest genetic association (Miesho et al. 2019). Two were on chromosome 8 (Vigun08g132300 and Vigun08g158000), but these genes differed from those found on chromosome 8 associated with early-exiting larvae in this study. Such differences among mapping studies are common (reviewed in Weiss 2008;Würschum 2012;Schielzeth et al. 2018), and might arise from differences in the genetic diversity captured in the mapping populations or in the specific ways that crop resistance was measured.
In mung bean (V. radiata), a close relative to cowpea, multiple recent studies have suggested that resistance to C. maculatus is mostly determined by a single major QTL (Chotechung et al. 2016;Kaewwongwal et al. 2017Kaewwongwal et al. , 2020. This QTL harbors genes for two polygalacturonase-inhibiting proteins (PGIPs) (Kaewwongwal et al. 2017;Zhang et al. 2021). Plant PGIPs are well known to be involved in defenses against pathogens; they inhibit polygalacturonases (PGs) that pathogens use to hydrolyze plant cell-wall pectins (De Lorenzo et al. 2001;Di Matteo et al. 2006). Recent evidence suggests that herbivorous insects also express PGs, and that PGIPs may protect plants from leaf-feeding insects as well (Haeger et al. 2020). Intriguingly, the PGIP-containing C. maculatus resistance QTL maps to 14.8-15.1 cM on chromosome 5 of the mung bean genome (Kaewwongwal et al. 2020) (for additional evidence that mung bean chromosome 5 underlies bruchid resistance, see Schafleitner et al. 2016). However, despite high overall synteny between mung bean and cowpea genomes (Lonardi et al. 2019a), we find no direct evidence of PGIPs on chromosome 5 in cowpea. Specifically, no PGIPs have been identified on chromosome 5, but PGIPs are found on other cowpea chromosomes. Nonetheless, one of our candidate genes, Vigun05g046000, is annotated as a "leucine-rich repeat protein", which is also true of PGIPs (Di Matteo et al. 2006). Evidence for a role of PGIPs in cowpea resistance might therefore emerge from refined annotations of the cowpea genome. At present, it appears that the genetic basis of resistance to C. maculatus differs between cowpea and mung bean, but further investigation of this conclusion is clearly warranted.
Selection for desirable crop traits would be simpler in the absence of epistasis because the beneficial effects of a particular gene variant would not strongly depend on genomic background. Nevertheless, few genomic studies have quantified the roles of epistasis vs. additive effects in accounting for variation in crop resistance to insects (Liu and Yan 2019;Soyk et al. 2020). On the whole, our results suggest that, for the cowpea genome, gene-gene interactions play at best a small role in explaining variation in the performance of seed beetles. We did obtain some statistical support for marginal epistasis among SNPs associated with development time, and these variants were in the same genomic region as those associated with additive-genetic variation. However, despite the evidence of marginal epistasis for these SNPs, including epistatic terms did not improve the explanatory power of the polygenic, genotype-phenotype model. This discrepancy raises a cautionary note for evaluating the importance of epistasis in plant resistance. In our analyses, adding epistasis terms actually caused a slight reduction in model's overall explanatory power. Given the hierarchical nature of the Bayesian models, adding epistasis terms may have slightly reduced the conditional probabilities of association for other individual loci. Evaluating the relative importance of epistasis is thus likely to require careful comparisons among multiple statistical models.

Implications for resistance assays and crop improvement
Comparisons between the second and third experiments revealed how intraspecific variation in pest life-histories can modify inferences about host resistance. For example, the CA and SI populations responded very similarly to genetic variation among RILs with respect to larval development time, but did not do so with respect to egg number. In a competitive environment, oviposition by SI females was likely strongly affected by the presence of eggs already on the seeds, and this factor may have outweighed or obscured variation among RILs in their intrinsic tendencies to elicit higher or lower amounts of egg-laying (Messina 1991). The number of eggs laid in either choice or no-choice tests is frequently measured as a component of grain-legume resistance to C. maculatus (e.g., Boeke et al. 2004;Cruz et al. 2016;Messina et al. 2018), but our results suggest that estimates of variation in host attractiveness can depend on the geographic population of the pest. Similarly, the number of eggs laid per 10 seeds was a good predictor of the number of emerging adults in the CA population, but contest competition within seeds caused a much weaker correlation in the SI population. Egg density and larval survival were negatively correlated as expected in the SI population, but not in the CA population.
Most importantly, the divergent behavior of SI and CA beetles ultimately led to a nearly fourfold difference in the total number of F1 adults (an overall measure of resistance), despite identical initial conditions of placing two pairs of beetles on each cowpea genotype for three days. Fortunately, our mapping analyses do not suggest that the two beetle populations use fundamentally different detoxification pathways, since there was a reasonable consistency in locations of causal variants for development time and survival. Future investigations of grain-legume resistance to C. maculatus may profit from using multiple pest populations and assay protocols. Cosmopolitan insect pests of stored products may commonly consist of genetically diverse "biotypes" as a consequence of founder effects, persistent genetic drift, or local natural selection (as is likely to occur if populations encounter different local hosts) (Downie 2010;Tuda et al. 2014;Semeao et al. 2012;Taggar and Arora 2017).
The cowpea MAGIC population has already been shown to contain wide phenotypic variation for diverse agronomic traits, including flowering time, growth habit, leaf shape and various seed characteristics, such as size and color (Huynh et al. 2018). Some of these traits provided evidence for strong transgressive segregation of the relevant alleles. For example, under conditions of drought stress, 11% of the RILs exhibited higher yield than any parent genotype. Significant transgressive segregation could therefore improve the likelihood of selecting lines that are agronomically superior to each of the elite parent cultivars (de los Reyes 2019). Across the three experiments in this study, a modest percentage of cowpea RILs similarly conferred slower development (13-20%), poorer survival to adult emergence (2-22%), and greater early-exit mortality ( ∼ 1-6%) than even the most resistant parental line. This study also identified reasonably high heritabilities and a relatively small number of causal variants for most key resistance traits, two factors that may facilitate combining major-effect alleles in breeding programs for crop improvement.
Further work is needed to examine genetic correlations between traits that appear to confer resistance to beetles in storage and those that mediate plant responses to abiotic and biotic stresses in the field (Lucas et al. 2012;Huynh et al. 2016). Interestingly, our resistance-associated region on chromosome 5 coincides with a known QTL that spans positions 3,748,359 to 4,572,228 in the cowpea genome and is associated with resistance to the cowpea aphid, Aphis craccivora (Huynh et al. 2015). These coincident QTLs suggest a possible pleiotropic gene that confers resistance to insects from very different ecological guilds: aphids feeding on phloem on whole plants vs. seed beetles feeding on cotyledons in dried seeds. It also remains to be seen whether seed traits that could significantly reduce the population growth of C. maculatus would be compatible with other important aspects of seed quality, which include size, shape, color, and texture (Huynh et al. 2018). For subsistence growers in regions such as sub-Saharan Africa, the development of cultivars that combine resistance to multiple stresses and market-preferred grain characteristics can reduce pesticide usage and yield a more reliable source of quality protein (Boukar et al. 2019;Horn and Shimelis 2020).