The UCR Minicore: a resource for cowpea research and breeding

Incorporation of new sources of genetic diversity into plant breeding programs is crucial for continuing to improve yield and quality, as well as tolerance to abiotic and biotic stresses. A minicore (the “University of California, Riverside (UCR) Minicore”) composed of 368 worldwide accessions of cultivated cowpea has been assembled, having been derived from the UCR cowpea collection. High‐density genotyping with 51,128 SNPs followed by principal component and genetic assignment analyses identified six subpopulations in the UCR Minicore, mainly differentiated by cultivar group and geographic origin. All six subpopulations were present to some extent in West African material, suggesting that West Africa is a center of diversity for cultivated cowpea. Additionally, population structure analyses supported two routes of introduction of cowpea into the U.S.: (1) from Spain to the southwest U.S. through Northern Mexico and (2) from Africa to the southeast U.S. via the Caribbean. Genome‐wide association studies (GWAS) narrowed several traits to regions containing strong candidate genes. For example, orthologs of the Arabidopsis FLOWERING LOCUS T lie within a major QTL for flowering time. In summary, this diverse, yet compact cowpea collection constitutes a suitable resource to identify loci controlling complex traits, consequently providing markers to assist with breeding to improve this crop of high relevance to global food and nutritional security.

It provides a major source of dietary protein, fiber, minerals, and vitamins for millions of people in sub-Saharan Africa (SSA), and fodder for livestock. Most of the production is in SSA by smallholder farmers, most of whom are women. It is also grown in many other parts of the world including Latin America, Southeast Asia, the Mediterranean Basin, and the United States (FAOSTAT, www.fao.org). Cowpea is well-known for its adaptation to heat and drought, and to soils with low fertility, making it a successful crop in arid and semi-arid regions where most other crops do not perform as well (Boukar et al., 2019).
However, breeding for increased heat and drought tolerance as well as for key agronomic traits and pest and disease-resistance is crucial as climate changes associated with global warming increase, and given that cowpea is primarily grown in regions that are quite vulnerable to climate change (Knox et al., 2012;Müller & Robertson, 2014).
The largest cultivar groups are probably Unguiculata and Melanophthalmus, which includes most grain and forage cowpeas, and Sesquipedalis, which is also known as "asparagus bean" or "yardlong bean" and is characterized by long and succulent pods that are consumed as a vegetable mostly in southeastern Asia (Boukar et al., 2019;Maréchal et al., 1978;. In addition to being grown for grain, cowpea is a source of nutritious fodder for livestock in dry savanna regions of sub-Saharan Africa (Dugje et al., 2009). One important aspect of cowpea agronomy is the time to flowering, as early-maturing types can in many cases be deployed as a strategy to capitalize on shortened periods of optimal growth, thus avoiding lateseason drought with its accompanying array of biotic stressors (Boukar et al., 2019). and developing core and minicore subsets to facilitate access to the diversity contained in the entire set of accessions is a common practice in most ex situ collections (Brown, 1995;Byrne et al., 2018).
Genetic and phenotypic evaluations of diverse collections are needed to fully utilize their potential in breeding programs. Several previous studies have reported on the genetic diversity of cultivated cowpea germplasm. Huynh et al. (2013) genotyped 442 cowpea landraces, predominantly from Africa, with a 1536-SNP GoldenGate assay (Muchero et al., 2009) and identified two major genepools, nominally West versus Southeast Africa. Later, 768 accessions mostly from the USDA cowpea collection were characterized using 5828 GBS (genotyping-by-sequencing) SNPs, revealing the presence of three major subpopulations (Xiong et al., 2016), again one from West Africa, with the other two subpopulations separating Asian, European, and some US accessions from those originating in India, Oceania, other parts of Africa, and the Americas. More recently, the majority of an IITA minicore collection (298 accessions) was genotyped using GBS with 2276 SNPs, also identifying three major subpopulations , but with more dispersion of West and Central African accessions across the three sub-populations. Carvalho et al. (2017) genotyped a smaller set of 96 worldwide cowpea accessions emphasizing Iberian Peninsula germplasm, but at a much higher SNP density using the Illumina Cowpea iSelect Consortium Array containing 51,128 SNPs . In this latter work, logical relationships were noted between the genetic composition of accessions and colonial-era movement of germplasm from the Iberian Peninsula to the Caribbean, and from Africa to South America. In general, it is evident that there is much yet to be clarified regarding the spread of cowpeas worldwide, and the extent to which certain genetic variants have taken hold in different regions at different times.
Here, we report the development of a minicore (henceforth referred as the "UCR Minicore") composed of 368 domesticated cowpeas selected from a larger set of $5000 accessions comprising the UC Riverside cowpea collection. This minicore was genotyped using the 51,128-SNPs Cowpea iSelect Consortium Array  and the genotypic information was used to gain a better understanding of the genetic diversity of domesticated cowpea. Although this is the first summary report and full SNP data release for the UCR Minicore, this material has been utilized for more focused work on seed coat color , seed coat pattern (Herniter et al., 2019), seed size , bruchid resistance (Miesho et al., 2019), plant herbivore resistance (Steinbrenner et al., 2020) and pod shattering (Lo et al., under review). This study evaluated additional traits of agronomic importance including flowering time, dry pod weight, dry fodder weight, and pod load score.
Genome-wide association studies (GWAS) were conducted, identifying significant SNPs and candidate genes associated with each of these traits.

| Plant materials and phenotyping
A total of 668 accessions representing a core subset from the University of California Riverside (UCR) cowpea germplasm collection were genotyped with an Illumina GoldenGate assay (Muchero et al., 2009) (1536 in previous studies (Huynh et al., 2013;Muchero et al., 2013). This SNP information, coupled with available phenotypic and passport data, was used to choose a smaller subset of nonredundant accessions, each of which was highly homozygous, collectively representing the genetic and phenotypic diversity of the larger core collection. A few additional accessions nominated by cowpea researchers and breeders based on regions of origin not represented among the core collection and possessing some specific traits were also included. The minicore of 368 accessions includes landraces and breeding materials from 50 countries in Africa, Asia, North and South America, Europe, and Australia (Table S1). Individual plants were grown from each of the 368 accessions in a greenhouse at UCR for genotyping (see Section 2.2) and seed production. Subsequent seed increases for distribution and phenotypic evaluations used seeds directly descended from these genotyped plants. was scored as the number of days from planting to when 50% of plants had at least one flower opened. Pod load score was recorded at plant maturity using a 1-3 scale with 1 for high pod load (80-100% of peduncles had 2-3 pods per peduncle), 2 for moderate pod load (80%-100% of peduncles had 1-2 pods per peduncle) and 3 for poor load score (80%-100% of peduncles had 1 or less pod per peduncle). At maturity, dry pods were harvested, and the other above-ground parts of the plant (leaves, twigs, and stems) were cut and rolled up. Both the pods and the fodder were sun dried for 1 week to determine dry pod weight and fodder weight respectively.

| SNP genotyping and data curation
Genomic DNA was extracted from young leaves of individual plants  , on the 368 minicore samples (Table S2).

| Genetic diversity, population structure and linkage disequilibrium analyses
Expected heterozygosity (He) and nucleotide diversity (π) values were calculated for all 48,425 SNPs in the minicore as a measure of genetic diversity. He was calculated asHe ¼ 1 À P k i¼1 Pi 2 , where P i is the frequency for the ith allele among a total of k alleles. π was evaluated as in .
The admixture model implemented in STRUCTURE v2.3.4 (Pritchard et al., 2000) was used to infer population structure of the UCR Minicore. Only SNPs with minor allele frequency (MAF) ≥ 0.05 were included in the analysis. STRUCTURE was first run three times for each hypothetical number of subpopulations (K) between 1 and 10, with a burn-in period of 10,000 followed by 10,000 Monte Carlo Markov Chain (MCMC) iterations. LnP(D) values were plotted and ΔK values were calculated according to Evanno et al. (2005) to estimate the optimum number of subpopulations. Plots were generated with Structure Harvester (Earl, 2012) ( Figure S1). Then, a new run using a burn-in period of 100,000 and 100,000 MCMC iterations was conducted for the estimated K to assign accessions to subpopulations based on a membership probability ≥0.80. In addition, principal component analysis (PCA) and linkage disequilibrium (LD) analysis were conducted on the same SNP set using TASSEL v5 (Bradbury et al., 2007).
Linkage disequilibrium (LD) was estimated for each chromosome as the correlation coefficient r 2 between pairs of SNPs, using SNPs with MAF ≥ 0.05 (42,711 SNPs). The decay of LD over physical distance was investigated by plotting pair-wise r 2 values and generating a locally weighted scatterplot smoothing (LOESS) curve in R. LD decay distance was determined when r 2 fell to the critical threshold estimated from the 95th percentile r 2 distribution for unlinked markers (r 2 = 0.15).

| GWAS
The mixed linear model (Zhang et al., 2010) implemented in TASSEL v.5 (Bradbury et al., 2007) (http://www.maizegenetics.net/tassel) was used for GWAS, with a population structure matrix (for K = 6) and a kinship coefficient matrix accounting for population structure and the relatedness of accessions, respectively. A false discovery rate (FDR; Benjamini & Hochberg, 1995) threshold of 0.01 was used to identify significant associations. The percentage contribution of each SNP to the total phenotypic variation was calculated using marker R 2 values from TASSEL multiplied by 100. Candidate genes within significant regions were identified from annotations of the cowpea reference genome v1.1 (https://phytozome.jgi.doe.gov/; Lonardi et al., 2019).

| Diversity and linkage disequilibrium in the UCR Minicore
The UCR Minicore was developed to represent available genetic and phenotypic diversity of cultivated cowpea while maintaining a sample size that can be managed by most researchers and breeders for evaluating traits of interest. This minicore is composed of 368 accessions from 50 countries including 242 landraces, 98 breeding lines, three accessions categorized as "weedy," and 25 accessions that are uncategorized (Table S1). It largely overlaps with the IITA minicore, with 233 accessions in common (Table S1), at least by name though not necessarily by genetic identity. The UCR Minicore showed great phenotypic diversity in pod and seed types and colors, flowering time and maturity, leaf shape, and plant architecture, among other traits. Expected heterozygosity (He) and nucleotide diversity (π) were calculated and both averaged 0.313 (the maximum for a biallelic SNP is 0.5).
The LD decay distance was also investigated for each cowpea chromosome using SNPs with MAF ≥ 0.05. The decay of LD with physical distance differed among chromosomes and ranged from 809 kb on Vu05 to 4,705 kb on Vu10 ( Figure S2). In addition to Vu05, the LD decay distance was shortest in Vu08 (813 kb) and Vu07 (1048 kb), while Vu01 and Vu04 had the largest decay distances after Vu10 (3803 and 3286 kb, respectively; Figure S2). On average there was one SNP per 14.9 kb, indicating sufficient marker density for GWAS.

| Population structure of the UCR Minicore
A total of 42,711 SNPs with MAF ≥ 0.05 were used to evaluate population structure in the UCR Minicore. STRUCTURE (Pritchard et al., 2000) was run for K = 1-10 and the inspection of both the estimated log probabilities of the data and the ΔK values calculated as in Evanno et al. (2005) supported the presence of six genetic subpopulations ( Figure S1). Accessions were assigned to each subpopulation using a membership coefficient ≥0.8 (accessions with membership coefficients less than 0.8 were considered "admixed"; Table S1). PCA F I G U R E 1 Phenotypic diversity in the UCR Minicore. Examples of diversity in (a) pod color and morphology and (b) seed coat color and pattern also showed a clear separation of the six subpopulations on the first three principal components ( Figure S3).
Subpopulation 1 included 31 accessions mostly from West Africa, 29 of which are landraces (Table S1) (Table S1). The largest subpopulation was Subpopulation 6 (69 accessions), which was composed primarily of landraces from Southeastern Africa. The remainder of the accessions (165) were considered "admixed" (Table S1). to the lack of passport information about the state of origin for most of them. Their cultivar names, however, traced to U.S. southern states (Table S1). The figure graphically shows that while all subpopulations are present in West Africa, three of them (1, 2, and 5) are predominant. Accessions from Southeastern Africa have ancestry primarily from Subpopulation 6, while most germplasm from Mediterranean countries belonged to Subpopulation 3. Interestingly, a closer look at the germplasm from the USA indicates that accessions from California have ancestry mostly from Subpopulation 3, while accessions from other U.S. states are predominantly Subpopulation 6 ( Figure 2). Subpopulation 4 is primarily from countries in Asia and Oceania.

| GWAS of agronomic traits
GWAS was conducted for four different agronomic traits including days to flowering (DTF), pod load score, dry pod weight, and dry fodder weight using 42,711 SNPs with MAF ≥ 0.05 (see Section 2). DTF was evaluated in five different environments in the USA and Nigeria, three of which are considered short-day environments (<12 hr of daylight) while the other two are considered long-day environments (>12 hr of daylight) (Table 1 and Section 2). Pod load score, dry pod weight, and dry fodder weight were evaluated in one environment in F I G U R E 2 Population structure in the UCR Minicore. The geographical distribution of the accessions belonging to the minicore is shown, with pie charts representing the average proportions from each genetic subpopulation for samples in each country. Pie chart sizes are proportional to the number of samples in each country and range from 1 to 80 accessions. Accessions within the U.S. were further divided as "California" and "other states/unknown," as most accessions in this second group did not have state information but their cultivar names traced to US southern states. The plot of ancestry estimates for K = 6 is shown at the bottom, with each bar representing the estimated membership coefficients for each accession in each of the six subpopulations (represented by different colors) T A B L E 1 Peak SNPs associated with DTF under five different environments. Asterisk on significant region (bp) indicates previously known QTL for DTF, as discussed in the text Nigeria. Significant marker-trait associations were identified for all traits and environments (Figures 3 and 4 and Tables 1 and 2 and S3 and S4).

| Days to flowering
Considering all five environments, 91 significant SNPs at 40 genomic regions were identified for DTF ( Figure 3 and  (Figure 3 and Tables 1 and S3). The phenotypic variation of significant SNPs ranged from 5% to 9% (Tables 1 and S3).
Five of the significant genomic regions corresponded to previously reported QTLs for DTF ( Figure 3 and  higher-density genotyping than previously, together with a larger and more diverse set of accessions.
Several candidate genes stood out for these five QTLs, as follows.
One of the previously reported QTL was identified in Thermal 2016 and is located at the beginning of Vu05 (371,674 bp;  found. This gene is orthologous to the Arabidopsis BRASSINAZOLE-RESISTANT 1 (BZR1), which positively regulates the brassinosteroid signaling pathway (He et al., 2002). The third QTL of interest was also identified in Riverside in 2017; it was located on Vu03 between 9,053,619 and 9,060,012 bp (Table 1 and Figure 3). The cowpea gene Vigun03g104200, which is an ortholog of the soybean E6, a main flowering and maturity gene, was identified 21-kb upstream of the peak SNP (2_24857).  Figure 4, and Table S4). Interestingly, the major QTL for dry fodder weight coincides with the QTLs for pod load score and dry pod weight (Table 2; Figure 4; Table S4). The colocation of these QTLs together with the correlations between the three traits suggests pleiotropic effects of a single gene or the existence of closely linked genes.

|
Genes within this common QTL region which spans 125 kb were

| Population structure and historical crop dispersal
The UCR Minicore and its associated high-density SNP data (51,128 SNPs) constitute a powerful combination of material and information resources to support genetic diversity analyses of cultivated cowpea.
Six subpopulations were identified in this minicore, all of which were represented to at least some extent in West African material ( Figure 2). This indicates that West Africa is a center of diversity for cultivated cowpea, as suggested by previous studies Padulosi & Ng, 1997;Steele, 1976). Five of the six subpopulations were composed mostly of landraces, while Subpopulation 2 included cowpea lines only from IITA's cowpea breeding program.
Based on how the six subpopulations split at different K numbers, as  (Table S1). Other work reported that Subpopulation 1 had much more pod shattering than Subpopulation 5 accessions (Lo et al., under review). In addition, all accessions in Subpopulation 1 have smooth seed coats, while most in Subpopulation 5 have rough coats (data not shown). Seed coat texture is an important quality trait that influences seed end-use; rough seedcoated varieties are preferred for food preparations requiring seed coat removal (e.g., making of Akara) as rough seed coats can be easily F I G U R E 4 Manhattan plots of genome-wide association studies (GWAS) on pod load score, dry pod weight, and dry fodder weight. Àlog10 (p values) are plotted against physical positions on the cowpea reference genome . The dashed red line in each plot indicates the 0.01 FDR-corrected threshold (4.62 for dry fodder weight and 4.63 for both pod load score and dry pod weight) T A B L E 2 Peak SNPs associated with dry fodder weight, dry pod weight, and pod load removed after soaking, while varieties with smooth seed coats are often preferred when cowpea is consumed as boiled intact seed (Singh & Ishiyaku, 2000). Also, rough seed coat types imbibe water quicker and generally have reduced cooking times compared to smooth seed coat types.  (Castetter & Bell, 1942;Herniter et al., 2020). In the Southeastern United States, cowpea seems to have been brought on slave ships, perhaps as provisions. The two distinct introduction routes of apparently genetically distinct cowpeas contrast with older studies predating genotyping, which have often assumed a single introduction.

Genotyping of the UCR
This conjecture of more than one route of cowpea introductions into the USA, and genetic distinctions between them that now are evident, is further supported by the findings of Carvalho et al. (2017), which showed relationships between a Cuban and Mediterranean accessions, and between sub-Saharan African and South American accessions, which represent yet another route of colonial-era dispersal of cowpeas.
Although the germplasm utilized in this study largely overlaps with that utilized by Huynh et al. (2013), only two genetic clusters were identified in that previous study which correspond to the two major "West Africa" and "Southeastern Africa" gene pools. As shown by Herniter et al. (2020) utilizing this same UCR Minicore material, at K = 2 Subpopulation 6 (Southeastern Africa gene pool) splits from the rest of the subpopulations, confirming that as a primary genetic differentiation between subpopulations of domesticated cowpeas. The smaller number of SNPs that were available for Huynh et al. (2013) to conduct population structure analyses, together with an emphasis on African landraces among accessions genotyped, most likely precluded further subdivision of the West African subpopulation. Carvalho et al. (2017) has been the only previous study to genotype a set of cowpea accessions at a high density (51,128 SNPs).
Although their focus was on genetic diversity of Iberian Peninsula cowpeas, results from that study largely agree with the findings reported in the present work. In particular, Subpopulations 1, 2, and 3 from the study of Carvalho et al. (2017), correspond to Subpopulations 4, 3, and 6 reported here, respectively. Their fourth subpopulation was composed only of four accessions from West Africa. This small representation of germplasm from West Africa certainly would preclude the detection of additional subpopulation structure present in that region.

| Flowering time
Flowering time is one of the most important agronomic traits, which affect environmental adaptation and yield potential. It is controlled by multiple genes via different pathways, and it is influenced by environmental conditions (Fornara et al., 2010).
Cowpea is generally a short-day plant and, although some accessions are day-neutral (photoperiod-insensitive), many cowpea genotypes are photoperiod sensitive and show a delay in flowering under long day conditions (Craufurd et al., 1996;Ehlers & Hall, 1996). The degree of photoperiod sensitivity can vary between accessions and is influenced by temperature (Ehlers & Hall, 1996).  (Fornara et al., 2010;Kobayashi et al., 1999 ;Yamaguchi et al., 2005). They share a similar mode of regulation, and overexpression of both FT and TSF results in precocious flowering (Kobayashi et al., 1999;Yamaguchi et al., 2005). The identification of multiple copies of FT genes in the cowpea reference genome  suggests that copy number variation could play an important role in the regulation of flowering time in cowpea, as reported for other crops (Díaz et al., 2012;Nitcher et al., 2013). Interestingly, in the cold-season legume chickpea a cluster of three FT genes has been identified under a major QTL controlling flowering time under shortday conditions (Ortega et al., 2019). Authors showed a collectively higher expression of the FT genes in the early-flowering domesticated lines respect to the late-flowering wild chickpea. Also, the flowering time QTL co-located with QTLs for growth habit and branching index, suggesting a possible involvement of FT genes in plant architecture as seen in other crops including Medicago truncatula (Laurie et al., 2011).
Another main flowering time locus identified under short days was located near Vigun04g057300, which is an ortholog of the Arabidopsis gene EID1. In Arabidopsis, EID1 mutations caused alterations in flowering induction (Marrocco et al., 2006). EID1 encodes an F-box protein that is a negative regulator in phytochrome A (phyA)specific light signaling (Dieterle et al., 2001). Mutations in EID1 causing the deceleration of the circadian clock have been selected during tomato domestication (Müller et al., 2016). Phytochrome modulation by environmental factors to control flowering time appears to be widespread in plants including cowpea (Mutters et al., 1989), where genetic variability related to EID1 potentially has been an important component of domestication.
A third QTL identified under a short-day environment on Vu05 was near Vigun05g077400, which is an ortholog of the Arabidopsis AGAMOUS-like 20 (AGL20) gene. AGL20 is an integrator of different pathways controlling flowering and is considered a central component for the induction of flowering (Borner et al., 2000;Lee et al., 2000).
Overexpression of AGL20 in Arabidopsis suppresses late flowering and delays phase transitions from the vegetative stages of plant development (Borner et al., 2000;Lee et al., 2000).
Under the long-day conditions of Riverside (CA, USA), most of the accessions in the UCR Minicore showed a delay in flowering, and about one fourth of the minicore did not flower (Table S1). Under long days, a major flowering time QTL was identified on Vu09 near Vigun09g059700, an ortholog of the Arabidopsis AGL8/FUL gene.
Vigun09g059700 is also an ortholog of the common bean gene Phvul.009G203400, which is a candidate for a photoperiod response QTL (DTF-9.5) identified in an Andean RIL population (Gonz alez et al., 2021). AGL8/FUL encodes a MADS-box family transcription factor that regulates inflorescence development and is negatively regulated by APETALA1 in Arabidopsis (Mandel & Yanofsky, 1995). In addition, AGL8/FUL promotes floral determination in response to farred-enriched light (Hempel et al., 1997). Loss of AGL8/FUL function in Arabidopsis caused a delay in flowering time (Ferr andiz et al., 2000;Melzer et al., 2008), suggesting that the cowpea ortholog of FUL (Vigun09g059700) could similarly be involved in the response to photoperiod.
Another region associated with flowering under long day was identified near genes Vigun11g169600 and Vigun11g169400. Both genes encode AP2/B3-like transcription factors that are orthologs to the Arabidopsis VRN1 gene. VRN1, a close homolog of FUL, promotes flowering after prolonged cold (vernalization) in different plant species (Levy et al., 2002;Lü et al., 2015;Putterill et al., 2004;Sung & Amasino, 2005;Yan et al., 2003). It is intriguing to identify a vernalization-pathway gene in a warm-season legume that does not need to undergo vernalization before flowering. However, vernalization pathway genes have been identified in other warm-season plants including soybean (Lü et al., 2015). In particular, Glyma11g13220, which was a homolog of the Arabidopsis VRN1, was responsive to photoperiod as well as to low temperatures in soybean. This vernalization pathway gene could also be functional in cowpea.
The cowpea gene Vigun09g244300, located on Vu09 and encoding a protein of the BES1/BZR1 family, is another strong candidate gene for days to flowering under long-day conditions.
Vigun09g244300 is an ortholog of the Arabidopsis BRASSINAZOLE-RESISTANT 1 (BZR1), one of the main regulators of the brassinosteroid signaling pathway (He et al., 2002). In Arabidopsis, brassinosteroid signaling inhibits the floral transition and promotes vegetative growth. Furthermore, brassinosteroid-deficient mutants cause a strong delay in days to flowering (Li et al., 2018). Lastly, another QTL was found on Vu03 near Vigun03g104200, an ortholog of the soybean E6 gene, which affects both flowering and maturity (Li et al., 2017;Sedivy et al., 2020). E6 plays an important role in the long-juvenile trait (delayed flowering) in soybean (Bonato & Vello, 1999).

| Plant productivity traits
Pod load score, dry pod weight, and dry fodder weight are important traits related to plant productivity. A cowpea genotype with high pod load (i.e., low pod load score) demonstrates high grain yield potential.
High pod load is a result of a high number of pods per plant, which is an indication of low rate of flower abortion. Generally, low flower abortion is associated with resistance to insect attack (typically flower thrips or maruca) and high night temperatures (>18 C) (Patel & Hall, 1990). Negative correlations have been identified between pod load score and dry pod weight, as well as between pod load score and number of pods per plant, and between pod load score and grain yield (Garcia-Oliveira et al., 2020). The negative correlation between pod load score and grain yield holds true as long as no attack by pod sucking bugs occurs, and in such a situation selection of high grain yielding genotypes can be made using pod load score. However, there are high positive correlations between dry pod weight and grain yield. In some studies, dry pod weight has been negatively correlated with dry fodder weight (Samireddypalle et al., 2017;Singh et al., 2003) except for some lines with dual purpose characteristics (Boukar et al., 2016;Timko & Singh, 2008).
We identified one major QTL associated with pod load score, dry pod weight, and dry fodder weight. Seven genes encoding members of the CNGC family proteins were located within the QTL region: four of them (Vigun04g039300, Vigun04g039400, Vigun04g039800, and Vigun04g039900) are orthologs of the Arabidopsis CNGC20 gene, also called CNBT1 (cyclic nucleotide-binding transporter 1), which is involved in the response to nematodes (Jha et al., 2016;Kaupp & Seifert, 2002). The other three (Vigun04g039500, Vigun04g039600, and Vigun04g039700) are orthologs of the Arabidopsis CNGC19 gene, which is involved in herbivore response (Jha et al., 2016;Meena et al., 2019). Interestingly, this region corresponds to the major Rk locus for root-knot nematode resistance identified in cowpea Ndeve et al., 2019). However, the authors are not aware of any nematode infestation or herbivore damage at the Minjibir field location. Members of the CNGC family of proteins have also been involved in plant tolerance to heavy metals (Moon et al., 2019).
CNGC20 and CNGC19 are Ca 2+ permeable channels, which play essential roles in the regulation of plant immunity and the response to abiotic stresses and thereby may influence plant productivity.
Interestingly, three CNGC calcium channels (a, b, and c) are involved in nodulation in Medicago truncatula (Charpentier et al., 2016). They are needed for inducing the oscillations in calcium concentrations that mediate plant responses to rhizobial bacteria (Roy et al., 2020). Therefore, a role of these genes in cowpea nodulation is also plausible. Further studies are necessary clarify any possible role of these cowpea homologs in regulating pod load score, dry pod weight, and dry fodder weight.

ACKNOWLEDGMENTS
The authors acknowledge Jasmine Dixon for assistance with retrieval

ETHICS STATEMENT
This study does not require any ethical approval.