Genetic Diversity and Genome-Wide Association in Cowpeas (Vigna unguiculata L. Walp)

Abstract

Cowpea is one of the most popular dry-land legumes cultivated for food and forage in arid and semi-arid areas. Genetic diversity for global germplasm can be organized into core collections providing optimum resources to serve breeding requirements. Here, we present diversity analysis and genome-wide association study (GWAS) results for part of the cowpea core collection of the United States Department of Agriculture (USDA) along with breeding line controls. Included in the analysis were a total of 373 accessions analyzed with 6880 Single Nucleotide Polymorphism (SNP) markers from Genotyping by Sequencing (GBS). Population structure differentiated accessions into two groups irrespective of geographical origin and formed three clusters based on taxa upon phylogenetic analysis. A total of 56 SNPs were significantly associated to nine traits including pod length (25 Quantitative Trait Nucleotides, QTNs), seed anti-oxidant content (7 QTNs), dry pod color (7 QTNs), plant maturity (5 QTNs), flower color (5 QTNs), seed weight (4 QTNs), tolerance to low phosphate (1 QTN), growth habit (1 QTN), and response to rock phosphate (1 QTN) using Bayesian-information, Linkage-disequilibrium Iteratively Nested Keyway (BLINK), and Fixed and random model Circulating Probability Unification (FarmCPU) association models. Key genes related to all significant SNPs were identified based on annotations of the cowpea reference genome, including a flavonoid gene controlling flower color (Vigun08g040200.1), a root nodulation regulator for tolerance to low phosphate (Vigun11g168000.1), and numerous genes involved in signaling, biosynthesis, metabolite transport, and abiotic stress. Our results highlight the importance of maintaining public phenotyping databases at USDA and strengthening collaborations for data collection in cowpea to maximize research impacts.

1. Introduction

Cowpea (Vigna unguiculata L. Walp) is an important legume cultivated globally for food, forage, vegetables, and soil nutrient enhancement. In addition to its nutritive value, cowpea is drought tolerant as a legume and one of the most prominent pulse crops in arid and semi-arid regions [1,2], including the Sahel and Savannah regions of Africa, where it is one of the only high-nitrogen-fixing legumes to improve soil fertility when in symbiosis with Rhizobia spp. bacteria [3,4]. Nigeria is the largest producer of cowpeas, followed by several other countries in Sub-Saharan Africa (FAOSTAT, 2023, https://www.fao.org/faostat/en/#data/QCL, accessed on 14 February 2024). Additional areas of traditional cowpea production include the Mediterranean, the Caribbean, and the Southeastern United States. Cowpeas have high agricultural potential in Eastern and South Asia and many drier parts of Latin America or North America as global warming worsens, making it one of the main legumes for climate-change adaptation.

Cowpeas have tolerance to various climate and soil-based abiotic stresses such as high temperatures and salinity [5]. Some biotic constraints are concerning, but generally, few bacterial and fungal diseases occur on cowpeas relative to most legumes [6]. They are plagued by Striga weeds, by viruses borne by insect vectors, and by many arthropod pests, especially those in storage, but tend to be robust crops with few field losses [7].

There is moderate genetic diversity in extant cowpea germplasm, with greater differentiation present between regions compared to within countries [8]. In our previous work, we found that vegetable cowpeas (known as yardlong beans) in Asia are very distinct from the grain and fodder cowpeas grown predominantly in Africa [9]. The development of diverse sets of representative core and mini-core collections in cowpea has enhanced delineation of population structure and genetic diversity while also informing germplasm conservation efforts and breeding strategies [10,11,12,13]. Similarly, availability of genetic tools, especially single nucleotide polymorphisms (SNPs) markers, has boosted the accuracy of linkage map construction and characterization of genetic mechanisms influencing important traits in cowpea through quantitative trait loci (QTL) mapping or genome-wide association studies (GWAS). For instance, the Illumina Cowpea iSelect Consortium Array has been widely used in numerous cowpea GWAS studies, as reviewed by various authors [14,15]. This includes genetic mapping in California and West Africa for seed size [16], resistance to Fusarium wilt [17], pod length [18], plant immunity to herbivores [19], drought tolerance [20], salt tolerance [21], aphid resistance [22], and flowering time [23].

On the other hand, multiple cowpea GWAS publications have used Genotyping by Sequencing (GBS) for de novo SNP discovery, as it is a germplasm genotyping platform [24] that varies based on the restriction enzymes involved, as reviewed by Blair et al. [14]. For example, the GBS method has been deployed by breeding groups in Arkansas and Texas for genetic mapping of various traits in cowpea, including resistance to bacterial blight [25], salt tolerance in seedlings [26], plant growth habit [27], aphid resistance [28], resistance to cowpea mosaic virus (CPMV) [29], and iron deficiency chlorosis [30], but all analyses were conducted without a full reference genome until recently.

Our interest in the current study was to use the GBS method, together with the newly released reference genome, to conduct a GWAS study to identify loci involved in novel traits evaluated on the USDA cowpea core collection. We present genomic regions associated with various seed, pod, and floral traits in a cowpea mini-core collection. These traits are key targets in cowpea breeding as they relate to yield (growth habit, pod length, and seed weight), nutrition (antioxidant content), tolerance to abiotic stress (early maturity, tolerance to low phosphate, and efficient use of rock phosphate), and phenotypic markers (flower and pod color).

Previous GWAS studies with some of these same traits and similar groups of USDA germplasm or new mapping populations have found significant SNPs. For example, plant growth habit in a diversity panel of 487 cowpea accessions detected 10 significant SNPs [27]. A total of 17 significant SNPs were associated with antioxidant content in cowpea using three statistical models, but only two were associated across models [31]. For pod length, 72 marker-trait associations (MTAs) were identified in GWAS [18]. However, no SNPs were associated with dry pod color [32]. Of the pod MTAs, 55 could be located on the cowpea genetic linkage map, but only two from chromosome Vu03 and Vu08 were associated with pod length in a recombinant inbred lines (RIL) population of a biparental cross between a domesticated (IT99K-573-1-1) and a wild (TVNu-1158) cowpea accession [33]. Three QTLs associated with 100 seed weight located on Vu01, Vu06, and Vu08 were detected from this same population [33] and two further MTAs for seed weight were identified in a recently developed Multi-parent Advanced Generation Inter-Cross (MAGIC) cowpea population [34], where two MTAs for plant growth habit on Vu01 and Vu09 were also found. In terms of abiotic stress tolerance, adaptation to low phosphorous conditions was associated with 10 SNPs in previous mapping [35].

As mentioned earlier, most of these GWAS reports occurred prior to the availability of a fully annotated cowpea genome and the physical locations of significant SNPs were unknown. This complicates the design of applicable SNP markers from these previous studies for marker-assisted selection of the corresponding traits in cowpea breeding programs. Furthermore, any functional genes underlying these genomic regions were also unknown and remain uncharacterized. Therefore, analysis with genome-confirmed GBS markers is merited.

The objectives of this work were to use new genotyping of the USDA core collection to (1) discover SNP markers derived from GBS with well-characterized physical genome placement information; (2) characterize genetic diversity within the 373 cowpea accessions from the USDA based on the new GBS markers; and (3) identify significant SNPs and corresponding quantitative trait nucleotides (QTNs) through GWAS analysis for nine important agronomic, plant physiological, and domestication traits in cowpea.

2. Materials and Methods

2.1. Plant Material

A total of 373 genotypes from three cowpea subspecies were used in this study (Supplementary Table S1) including 295 plant introductions (PIs) of grain cowpea (V. u. ssp. unguiculata), 26 accessions of yardlong bean (ssp. sesquipedalis), and 13 accessions of forage cowpea (ssp. cylindrical). Another 39 cowpea cultivars were parents from two cowpea breeding programs. The first group of grain cowpea, yardlong bean, and forge cowpea were provided by USDA Germplasm bank in Griffin, GA, while the parental cowpea lines were provided by University of California (UC), Riverside, and UC-Davis.

2.2. DNA Extraction and Genotyping by Sequencing (GBS)

Six seeds of each the 373 cowpea accessions were disinfected by diluted bleach for ten minutes and geminated in a petri dish maintained in the dark growth chamber for four days. DNA was extracted from the cotyledon and the shoot apex of the cowpea sprouts using DNeasy Plant DNA miniprep kits (Qiagen, Hilden, Germany) according to the manufacturer’s procedures. The DNA quality and concentrations were detected by gel electrophoresis and FLUOstar Omega (BMG LABTECH, Cary, NC, USA). A total of 1.5 μg of DNA from each sample was sent to the Institute of Biotechnology, Cornell University for sequencing.

Genomic DNA was digested with the ApeKI restriction enzyme. Barcode adaptor ligation, sample pooling, and amplification for sequencing library construction were performed according to protocols described by Elshire et al. [24]. Single-end sequencing of the 95-plex library was performed with the Genome Analyzer II next generation sequencing platform (Illumina Inc., San Diego, CA, USA).

2.3. Phenotyping of Plant Traits

Phenotyping data of this collection were extracted from the USDA-ARS Germplasm Resources Information Network (USDA-GRIN) that was available at the descriptor webpage of the United States National Plant Germplasm System with the crop filter ‘VIGNA’ (https://npgsweb.ars-grin.gov/gringlobal/descriptors, accessed on 14 February 2024). Phenotyping with field or greenhouse trials was conducted by the national curator as the main observer from the late 1980s to 2012. He mostly used Griffin GA as his main location but occasionally collected data from St. Croix, Virgin Islands, and Isabela, Puerto Rico winter nurseries as listed in GRIN. A total of nine traits were considered for GWAS. Of these, five were quantitative in nature and evaluated numerically, while four were categorical and defined based on descriptors from USDA-GRIN.

The categorical descriptors [36] included: (1) Plant maturity, being the number of days required for pods to mature after sowing date and categorized as early (0–52 days), normal (52–104 days), late (104–156 days), and very late (156–208 days); (2) Growth habit, classified as short/erect, semi-prostrate, or prostrate/long vined; (3) Flower color evaluated at full bloom and coded as white, lavender, purple, yellow, pink, blue, or mixed; and (4) Dry Pod Color, categorized as straw, green, purple, black, speckled, or brown, with some genotypes having variable pod colors and being removed from analysis.

Quantitative traits included the following: (5) Pod length, measured on mature pods from the apex to the connection at the peduncle and averaged per genotype; (6) Seed weight, recorded as weight in grams (g) of 100 seeds for each genotype; (7) Antioxidant Activity (AOA), obtained by chemical analysis carried out in College Station, Texas and measured in micrograms trolox equivalents/gdw (gdw = grams of dry weight), based on the method of Koleva et al. [37] where 2 g of seed from each accession were ground, methanol-extracted, and treated with 2,2-Diphenyl-1-picrylhydrazyl (DPPH).

Two other traits were based on phenotyping the whole plant, namely (8) the growth response to low phosphorus (P) soil conditions measured in a greenhouse trial conducted in the Texas AgriLife Research Station in Bushland, Texas, using pots filled with Betis sand with pH 4.7 in 1:1 water: soil mixture, and plant-available P level of 3 mg kg⁻¹ and the Melich III method [38] used previously for GWAS [35]. After 8 weeks, plant heights were measured, and shoots were harvested, dried, and weighed. Total biomass (shoot dry weight plus root dry weight) was calculated. Results were grouped into five categories from efficient to inefficient at P uptake. No phosphate treatment was considered control. Trait (9) was response to rock phosphate, where 300 mg P kg⁻¹ was added, and eight weeks afterwards, shoots and roots were harvested, dried, and weighted; total biomass was again calculated and the PIs were categorized from highly responsive to low responsive.

2.4. SNP Identification

The raw sequence data were analyzed with the GBS discovery pipeline in TASSEL GBS V2 [39]. The FASTQ raw files and sample key files, with information of plate layout and bar codes for each genotype, were used to construct a GBS database for the identification of SNPs. Only the sequence reads containing barcode sequence followed by the sticky-end sequence of an ApeKI restriction enzyme cut site (GˇCWGC) were trimmed to 64 bases and stored in one database. Reads that did not match any barcode or cut site remnant were excluded from the analysis, as were reads containing unidentified bases (N) and reads with adapter dimers. Subsequently, the barcoded sequence reads with tags present more than three times were sorted and collapsed into unique sequence tags with position information, and then aligned with the Vigna unguiculata genome v1.1 [40] from Phytozome (https://phytozome.jgi.doe.gov/pz/portal.html, accessed on 14 February 2024) using the Burrows–Wheeler Aligner (BWA) algorithm [41]. All newly discovered SNPs were scored for coverage, depth (minimum 5×), genotypic information (minimum 80% at the loci and the genotype level), heterozygosity (maximum of 70%), and Minimum Allele Frequency (MAF of 0.05). The quality score of 10 was applied for the validation of any given locus. The resultant filtered marker set included 6880 SNPs and eliminated 17 genotypes with missing data, resulting in a matrix for 356 genotypes out of the 373 accessions used initially for DNA.

2.5. Population Structure and Genetic Diversity

Population structure was evaluated based on the allele calling using a Bayesian-based unsupervised clustering approach implemented in the software STRUCTURE 2.3.4 [42]. The K (population number) value was set from 1 to 10 with five replicates for each K value. Burn in and replicate number for each run were set up at 25,000 and 100,000, respectively, and all analyses were run through StrAuto v1.0 [43], a multi-threading python script. The best-fit K value was determined by generating a plot of the mean likelihood values per K using Evanno transformation in STRUCTURE HARVESTER v1.0 [43,44,45]. The distribution of each individual in each cluster was determined by CLUMPP v1.0 [46], and geographical origin was used to define the predominant cluster source for display in the dendrogram. Cluster analysis by Neighbor-Joining algorithm and principal component analysis (PCA) method [47] were performed to evaluate population structure and to clarify the spatial genetic relationships between accessions. PCA scores were displayed with DARwin v6.0 [48]. Analysis of molecular variance (AMOVA) was performed with polymorphic SNP markers using Arlequin v3.5 [49].

2.6. Population Stratification and Kinship

Random and fixed effects were estimated using the filtered 6880 SNP markers and 356 cowpea accessions to reduce the false positive rate of each GWAS model. Random effects comprised kinship relationships, while fixed effects accounted for population clustering. The kinship matrix was built using the VanRaden algorithm from the GAPIT package in R’s v. 4.1.2 environment [50]. On the other hand, the unsupervised population clustering was explored using the molecular principal component analysis from the same GAPIT v 3.0 software and, for each case, the first two and three components were plotted in R [51].

2.7. Genome-Wide Association (GWAS) Mapping

The 6880 SNP markers were inspected for associations with all nine traits surveyed across the 356 cowpea accessions remaining in the dataset after GBS marker evaluation. Associations were determined using the GWAS algorithms FarmCPU and BLINK as implemented in R’s GAPIT v3.0 package [51]. These models are known to boost GWAS’s statistical power while efficiently controlling the false-positive rate [52]. Kinship and each population stratification scenario (from K = 3 to K = 5) were respectively considered as random and fixed effects for a total of 24 models (four traits inspected with FarmCPU and BLINK models, considering, in each case, three population strata). Significant associations were inferred using a strict Bonferroni correction of p-value at an α = 0.05, implying an effective significance threshold of −log₁₀ (7.27 × 10⁻⁶) = 5.14 for all GWAS models. The rate of false positives was evaluated by visual inspection of the Q–Q plots. Circular Manhattan plots and their respective Q–Q plots were generated in RIdeogram [53].

3. Results

3.1. SNP Diversity, Population Structure, and Genetic Diversity of Cowpea Panel

Overall, 76,114 polymorphic SNPs were discovered across the 373 accessions/cultivars, with only 17 accessions with high missing data values (>70% threshold) excluded, resulting in a database for 356 genotypes. The data were subsequently filtered by removing SNPs with rare alleles (<5%), a high-missing ratio (>5%), or high heterozygosity (>70%), leading to a final dataset of 6880 markers for further analyses.

Population structure according to peak delta K separated all accessions into two groups (K = 2), irrespective of taxonomy or geographical origin (Figure 1a,b). The first group (Q1) consisted of 84 accessions (23.6% of total accessions), while the second group (Q2) had 117 accessions (32.9%) (Figure 1b). The remaining 155 accessions (43.5%) had varied levels of admixture from both Q1 and Q2. Accessions in Q1 comprised 79 landraces of three taxa (V. unguiculata, V. unguiculata ssp. unguiculata, and V. unguiculata ssp. cylindrica) from 26 countries and five breeding lines. Accessions in Q2 consisted of 14 breeding lines while the remaining 103 accessions belong to five taxa (V. unguiculata, V. unguiculata ssp. unguiculata, V. unguiculata ssp. sesquipedalis, V. unguiculata ssp. cylindrica, and V. unguiculata ssp. pubescens) from 30 countries. At K = 3, V. unguiculata ssp. sesquipedalis accessions were clearly separated from the other subspecies. Diversity analysis classified 356 accessions into three clusters, loosely based on their geographic origin.

Accessions from East and West Asia formed cluster I. The second cluster (cluster II) had two sub-clusters. Latin America, Europe, and Middle East comprised the first sub-cluster, while Central and West Africa accessions formed the second sub-cluster. Others from South and East Africa formed Cluster III. The scattered distribution across clusters for North American accessions showed their complex genetic background.

3.2. Genetic Diversity among Yardlong Types, Grain Cowpeas, and Their Wild Relatives

Although Neighbor-Joining analysis did not cluster the 356 accessions in three clusters based on their taxa, V. unguiculata ssp. sesquipedalis (yardlong bean), an important vegetable cowpea consumed in South and Southeast Asia, was clustered into Cluster I, being differentiated from accessions that belonged to other taxa from other countries. V. unguiculata ssp. cylindrica was distributed into both Cluster I and Cluster II. One accession that belonged to V. unguiculata ssp. pubescens from Tanzania in Cluster I, which has been used as a parent in cowpea breeding program, was found to be closely related to an accession that belonged to V. unguiculata ssp. unguiculata. The principal component analysis (PCA) also supported the separation of this set of accessions into three clusters, where yardlong bean accessions were closely grouped together while wild cowpea relatives (V. unguiculata ssp. cylindrica were clustered with grain cowpeas.

The variation within cowpea subspecies was supported by AMOVA analysis (

Table 1. Analysis of molecular variance based on three cowpea sub-species groups.

Source of Variation	Degrees of Freedom	Sum of Squares	Variance Components	Percentage of Variation
Variation partition (Among three sub-species groups)
Among groups	2	18,782.60	160.48	12.75
Among individuals within groups	75	126,940.08	594.23	47.21
Within individuals	78	39,318.50	504.08	40.05
Total	155	185,041.17	1258.79
Variation partition (Among Vigna unguiculata ssp. unguiculata and V. unguiculata ssp. sesquipedalis)
Among groups	1	15,703.60	219.09	15.36
Among individuals within groups	63	117,856.43	663.40	46.51
Within individuals	65	35,356.50	543.95	38.13
Total	129	168,916.55	1426.43
Variation partition (Among Vigna unguiculata ssp. unguiculata and V. unguiculata ssp. cylindrica)
Among groups	1	4654.56	72.87	5.94
Among individuals within groups	49	89,727.87	677.62	55.25
Within individuals	51	24,273.50	475.95	38.81
Total	101	118,655.93	1226.44
Variation partition (Among Vigna unguiculata ssp. cylindrica and V. unguiculata ssp. sesquipedalis)
Among groups	1	7524.12	174.74	15.52
Among individuals within groups	38	52,841.96	439.11	38.99
Within individuals	40	20,494.00	512.35	45.49
Total	79	80,860.08	1126.21

), which indicated that most of the variance occurred among groups, in other words, among subspecies, and accounted for 83.7% of the total variation, whereas 9.1% and 7.2% of the variation was attributed to differences within individuals and between sub-populations, respectively.

3.3. Genome-Wide Association Studies

A total of 56 significant QTNs were identified for all nine traits evaluated with both statistical models (

Table 2. Single nucleotide polymorphisms (SNPs) significantly associated with nine traits evaluated in the cowpea mini-core collection. An asterisk (*) indicates SNPs identified by both BLINK and FarmCPU models.

Trait	SNP	Chr	Position	p-Value	FDR-Adjusted p-Values	Model	SNP Effect
Antioxidant	CONTIG_3_347132	n/a	347,132	9.2 × 10−6	2.1 × 10−2	FarmCPU	176.81
Antioxidant	SVU02_25844432	VU02	25,844,432	2.8 × 10−9	9.6 × 10−6	BLINK
Antioxidant	SVU04_24674908	VU04	24,674,908	1.3 × 10−5	2.3 × 10−2	FarmCPU	101.15
Antioxidant	SVU04_39821223	VU04	39,821,223	5.0 × 10−7	1.7 × 10−3	FarmCPU	−142.33
Antioxidant	SVU08_36638094 *	VU08	36,638,094	2.8 × 10−8	6.3 × 10−5	BLINK
Antioxidant	SVU08_36638094 *	VU08	36,638,094	1.8 × 10−12	1.2 × 10−8	FarmCPU	−257.05
Antioxidant	SVU09_41120416	VU09	41,120,416	8.5 × 10−6	1.5 × 10−2	BLINK
Antioxidant	SVU10_20951283	VU10	20,951,283	2.6 × 10−9	9.6 × 10−6	BLINK
Dry pod color	SVU03_61354189 *	VU03	61,354,189	3.8 × 10−10	1.3 × 10−6	BLINK
Dry pod color	SVU03_61354189 *	VU03	61,354,189	1.5 × 10−6	4.2 × 10−3	FarmCPU	−0.46
Dry pod color	SVU03_9536356 *	VU03	9,536,356	8.4 × 10−7	9.5 × 10−4	BLINK
Dry pod color	SVU03_9536356 *	VU03	9,536,356	4.1 × 10−6	7.0 × 10−3	FarmCPU	−0.37
Dry pod color	SVU05_6973585 *	VU05	6,973,585	2.3 × 10−13	1.6 × 10−9	BLINK
Dry pod color	SVU05_6973585 *	VU05	6,973,585	4.9 × 10−8	3.4 × 10−4	FarmCPU	0.74
Dry pod color	SVU06_33534718 *	VU06	33,534,718	1.1 × 10−9	2.4 × 10−6	BLINK
Dry pod color	SVU06_33534718 *	VU06	33,534,718	1.9 × 10−6	4.2 × 10−3	FarmCPU	0.50
Dry pod color	SVU08_36154680	VU08	36,154,680	3.2 × 10−7	4.3 × 10−4	BLINK
Dry pod color	SVU09_16383249	VU09	16,383,249	8.0 × 10−6	1.1 × 10−2	FarmCPU	0.80
Dry pod color	SVU09_25649823	VU09	25,649,823	2.0 × 10−9	3.4 × 10−6	BLINK
Flower color	CONTIG_3_248443	n/a	248,443	5.5 × 10−6	9.0 × 10−3	BLINK
Flower color	SVU01_35872139	VU01	35,872,139	6.6 × 10−6	9.0 × 10−3	BLINK
Flower color	SVU03_31090324	VU03	31,090,324	1.0 × 10−6	3.5 × 10−3	BLINK
Flower color	SVU08_3874640	VU08	3,874,640	4.5 × 10−6	9.0 × 10−3	BLINK
Flower color	SVU09_5475451	VU09	5,475,451	7.5 × 10−7	3.5 × 10−3	BLINK
Growth habit	SVU04_32873146	VU04	32,873,146	2.8 × 10−8	1.9 × 10−4	BLINK
Maturity	SVU03_47710713	VU03	47,710,713	6.3 × 10−6	2.1 × 10−2	BLINK
Maturity	SVU04_7507872	VU04	7,507,872	3.3 × 10−6	2.1 × 10−2	BLINK
Maturity	SVU05_16327515	VU05	16,327,515	3.5 × 10−7	8.1 × 10−4	FarmCPU	8.96
Maturity	SVU09_12400987	VU09	12,400,987	1.2 × 10−9	8.5 × 10−6	FarmCPU	−11.55
Maturity	SVU09_16383398	VU09	16,383,398	3.5 × 10−7	8.1 × 10−4	FarmCPU	8.61
Pod length	CONTIG_3_141951	n/a	141,951	1.2 × 10−6	6.7 × 10−4	FarmCPU	4.74
Pod length	CONTIG_3_495433	n/a	495,433	9.6 × 10−8	4.1 × 10−5	BLINK
Pod length	SVU01_3813516 *	VU01	3,813,516	9.4 × 10−19	6.4 × 10−15	BLINK
Pod length	SVU01_3813516 *	VU01	3,813,516	5.3 × 10−14	3.6 × 10−10	FarmCPU	−5.62
Pod length	SVU01_3814443 *	VU01	3,814,443	8.3 × 10−8	3.8 × 10−5	BLINK
Pod length	SVU01_3814443 *	VU01	3,814,443	5.1 × 10−11	5.8 × 10−8	FarmCPU	−28.15
Pod length	SVU01_3814444 *	VU01	3,814,444	8.3 × 10−8	3.8 × 10−5	BLINK
Pod length	SVU01_3814444 *	VU01	3,814,444	5.1 × 10−11	5.8 × 10−8	FarmCPU	−28.15
Pod length	SVU01_3814447 *	VU01	3,814,447	8.3 × 10−8	3.8 × 10−5	BLINK
Pod length	SVU01_3814447 *	VU01	3,814,447	5.1 × 10−11	5.8 × 10−8	FarmCPU	−28.15
Pod length	SVU01_3814467 *	VU01	3,814,467	8.3 × 10−8	3.8 × 10−5	BLINK
Pod length	SVU01_3814467 *	VU01	3,814,467	5.1 × 10−11	5.8 × 10−8	FarmCPU	−28.15
Pod length	SVU02_22876112	VU02	22,876,112	8.9 × 10−9	7.6 × 10−6	BLINK
Pod length	SVU03_37127830	VU03	37,127,830	1.9 × 10−6	1.0 × 10−3	FarmCPU	−6.59
Pod length	SVU05_40768813 *	VU05	40,768,813	1.2 × 10−8	8.8 × 10−6	BLINK
Pod length	SVU05_40768813 *	VU05	40,768,813	2.0 × 10−8	1.5 × 10−5	FarmCPU	3.95
Pod length	SVU06_26010349	VU06	26,010,349	7.8 × 10−16	2.7 × 10−12	BLINK
Pod length	SVU06_563046	VU06	563,046	2.0 × 10−10	1.7 × 10−7	FarmCPU	−25.52
Pod length	SVU06_563048	VU06	563,048	2.0 × 10−10	1.7 × 10−7	FarmCPU	−25.52
Pod length	SVU07_32970747	VU07	32,970,747	1.2 × 10−13	1.7 × 10−10	BLINK
Pod length	SVU07_39887364	VU07	39,887,364	1.6 × 10−12	5.5 × 10−9	FarmCPU	−5.80
Pod length	SVU07_8419865	VU07	8,419,865	6.6 × 10−8	4.5 × 10−5	FarmCPU	3.47
Pod length	SVU08_2723206	VU08	2,723,206	4.5 × 10−6	2.0 × 10−3	FarmCPU	−3.73
Pod length	SVU09_29209316	VU09	29,209,316	1.8 × 10−9	1.7 × 10−6	BLINK
Pod length	SVU10_1257951	VU10	1,257,951	2.3 × 10−7	1.4 × 10−4	FarmCPU	4.29
Pod length	SVU10_32475355	VU10	32,475,355	9.0 × 10−11	1.0 × 10−7	BLINK
Pod length	SVU10_32507465	VU10	32,507,465	1.8 × 10−8	1.2 × 10−5	BLINK
Pod length	SVU10_33225256	VU10	33,225,256	2.5 × 10−6	1.2 × 10−3	FarmCPU	3.97
Pod length	SVU10_33358174	VU10	33,358,174	1.2 × 10−13	1.7 × 10−10	BLINK
Pod length	SVU11_36410645	VU11	36,410,645	3.8 × 10−8	2.4 × 10−5	BLINK
Pod length	SVU11_971859	VU11	971,859	6.8 × 10−14	1.6 × 10−10	BLINK
ResponseLow P	SVU10_3706971	VU10	3,706,971	9.9 × 10−8	6.8 × 10−4	BLINK
Rock phosphate	SVU11_37561353	VU11	37,561,353	2.7 × 10−6	1.9 × 10−2	BLINK
Seed weight	SVU03_46623413	VU03	46,623,413	2.6 × 10−7	8.9 × 10−4	BLINK
Seed weight	SVU03_50505045	VU03	50,505,045	7.7 × 10−6	1.3 × 10−2	BLINK
Seed weight	SVU03_52946137	VU03	52,946,137	2.9 × 10−10	2.0 × 10−6	BLINK
Seed weight	SVU04_16299944	VU04	16,299,944	1.6 × 10−6	3.6 × 10−3	BLINK

, Figure 2). Genes related to the QTNs were also identified and annotated. The most significant QTNs (n = 25) were associated with pod length, followed by dry pod color and antioxidant content (n = 7), and flower color and maturity (n = 5). Seed weight and maturity had four significant SNPs each, but only a single significant SNP was discovered for tolerance to low phosphate and rock phosphate. Overall, analyses using the BLINK model revealed significant SNPs for all nine traits while FarmCPU only identified significant SNPs in antioxidant content, dry pod color, maturity, and pod length. The number of SNPs identified by both methods for each trait was highly correlated (R² = 0.986).

3.4. Antioxidant Level

Of the SNPs associated with antioxidant content, SVU08_36638094 was significant in both the BLINK and FarmCPU models. SVU02_25844432, SVU09_41120416, and SVU10_20951283 were identified through BLINK while SCONTIG_3_347132, SVU04_24674908, and SVU04_39821223 were identified through FarmCPU. SVU08_36638094 is in the exon of gene Vigun08g202200 encoding a o-glycosyl hydrolases family 17 protein. SVU04_24674908 and SVU04_39821223 lie within Vigun04g104300, a leucine rich protein, and Vigun04g173800.1, an unknown protein. SVU02_25844432 is in a formin-like protein annotated as a tensin phosphatase C2 domain encoded by the gene Vigun02g103800. SVU09_41120416 is in the exon of Vigun09g241600, a transcript of unknown protein. SVU10_20951283 and SCONTIG_3_347132 were not located on any genes.

3.5. Dry Pod and Flower Color

Significant SNPs associated with dry pod color that were identified in both methods included SVU03_61354189, SVU03_9536356, SVU05_6973585, and SVU06_33534718. FarmCPU detected SVU09_16383249 on Chromosome 9 while SVU08_36154680 and SVU09_25649823 were detected by BLINK on chromosomes 8 and 9, respectively. SVU03_61354189 is located within Vigun03g406500, an ATPase domain, while SVU03_9536356 and SVU05_6973585 are located on zinc finger DOF binding transcripts Vigun03g106900 and Vigun05g076200, respectively. SVU06_33534718 was found in the exon of Vigun06g227800, an uncharacterized membrane protein.

All five significant SNPs obtained for flower color were detected through the BLINK model. SCONTIG_3_248443 fell within the exon of VigunL033300, a photosystem II P680 reaction center D1 protein (psbA). SVU01_35872139 was within the exon of an uncharacterized protein of gene Vigun01g177600. SVU03_31090324 and SVU08_3874640 were in an RNA polymerase, alpha chain C terminal domain and Vigun08g039500, an IQ calmodulin-binding motif, respectively. Gene Vigun09g054600 contains SVU09_5475451 and codes for a transporter endomembrane family protein 70.

3.6. Growth Habit and Maturity

One SNP significantly associated with cowpea growth habit was SVU04_32873146, located on Chromosome 4 within intron of the Vigun04g131400, a TUBBY-LIKE F-BOX PROTEIN 1-RELATED protein. In contrast, several significant QTNs were associated with maturity found by the BLINK model (SVU03_47710713 and SVU04_7507872). Each of these SNPs were located within introns of Vigun03g292300 that encode AT hook motif DNA-binding family protein and Vigun04g066800 for a sulphate adenylyltransferase catalytic domain protein. Three SNPs, SVU05_16327515, SVU09_12400987, and SVU09_16383398 were found to be significantly associated with the FarmCPU model and the former was located within the exon of Vigun05g136200, a Zinc-binding dehydrogenase. For the two SNPs on Chromosome 9, SVU09_12400987 fell within intron of Vigun09g089800, an NAD(P)-BINDING ROSSMANN-FOLD SUPERFAMILY PROTEIN, and SVU09_16383398 is in the exon of Vigun09g100100, a Photosystem II protein that transfers electrons within the cyclic electron transport pathway of photosynthesis activity.

3.7. Pod Length

Pod length had the highest number of significant SNPs (n = 25) detected for any of the traits in this study. The two analysis models uniquely detected 10 (BLINK) and nine (FarmCPU) significant SNPs, respectively, and a combined total of six SNPs similar to both models included five on Chromosome 1 (SVU01_3813516, SVU01_3814443, SVU01_3814444, SVU01_3814447, and SVU01_3814467) and one on Chromosome 5 (SVU05_40768813). SNPs SVU01_3813516, SVU01_3814443, SVU01_3814444, SVU01_3814447, and SVU01_3814467 were all within the transcript of an uncharacterized protein encoded by Vigun01g030600. All the other SNPs were distributed across nine chromosomes, with two in Chromosome 3. Chromosome 2 had one SNP SVU02_22876112 in the exon of Vigun02g076300 that encodes an ALDEHYDE DEHYDROGENASE FAMILY 2 MEMBER B7. Significant SNPs on Chromosome 6 included SVU06_26010349 in the exon of Vigun06g133700, a DEAD-BOX ATP-DEPENDENT RNA HELICASE 40 and two SNPs (SVU06_563046 and SVU06_563048) within the intron of Vigun06g001200, an ACYL-MALONYL CONDENSING ENZYME-RELATED protein. Both SVU07_32970747 and SVU07_39887364 were on Chromosome 7 within the exon of Vigun07g208100, an ABC TRANSPORTER G FAMILY MEMBER 25 and intron of Vigun07g285900, a photosynthetic PSBP-LIKE PROTEIN 2, respectively. Chromosomes 8 and 9 had one SNP each, SVU08_2723206 and SVU09_29209316, respectively, located in the exons of Vigun08g029800, activating signal co-integrator complex subunit 2 (ASCC2) and Vigun09g133300, a Threonine-specific protein kinase with a Leucine Rich Repeat (LRR_1). Chromosome 10 had five significant SNPs including SVU10_1257951, SVU10_32475355, SVU10_32507465, SVU10_33225256, and SVU10_33358174 in the respective exons of Vigun10g012100; a Ca²⁺-independent phospholipase A2, Vigun10g118300; a TRANSCRIPTION FACTOR TCP9, Vigun10g118800; a solute carrier family 35 (adenosine 3′-phospho 5′-phosphosulfate transporter, member B3 (SLC35B3, PAPST2), Vigun10g124000; a protein with unknown function and Vigun10g125200; and a short-chain alcohol dehydrogenase. Two SNPs from Chromosome 11, SVU11_971859 and SVU11_36410645, were located within the intron of uncharacterized protein encoded by gene Vigun11g008600 and the exon of Vigun11g154700, a Dof domain and zinc finger (zf-Dof), respectively. Four QTNs, SCONTIG_3_495433, SVU05_40768813, SVU07_8419865, and SVU11_36410645, were in non-coding regions, while SCONTIG_3_141951 was in the first exon of VigunL028900, a photosystem II P680 reaction center D2 protein (psbD).

3.8. Response Low P

The only significant SNP detected for tolerance to low phosphate was SVU10_3706971 and was within the first exon of Vigun10g029400 gene, a small subunit ribosomal protein S6e (RP-S6e, RPS6).

3.9. Rock Phosphate Response

Only one SNP (SVU11_37561353) was significantly associated with the response to rock phosphate. This was found in the exon of Vigun11g169100 encoding an Arginine decarboxylase/L-arginine carboxy-lyase protein.

3.10. Seed Weight

All four SNPs with significant association to seed weight were identified in the BLINK model and included three on Chromosome 3 (SVU03_46623413, SVU03_50505045, and SVU03_52946137) and one on Chromosome 4 (SVU04_16299944). The SNP SVU03_46623413 fell in a non-coding region, but SNPs SVU03_50505045 and SVU03_52946137 were in the intron of Vigun03g312100, a phosphatidylinositol glycan, class U (PIGU), and the exon of Vigun03g331800, a UBIQUITIN SPECIFIC PROTEINASE, respectively.

3.11. QTN Co-Localizations

Co-localizations of significant SNPs for various traits were observed, suggesting pleiotropic genetic regions. For instance, SVU03_46623413 (seed weight), SVU03_31090324 (flower color), SVU03_47710713 (maturity), and SVU03_37127830 (pod length) are within 4 Mbp of each other on Chromosome 3. On Chromosome 4, SVU04_16299944 (seed weight), SVU04_24674908/SVU04_39821223 (antioxidant content), SVU04_32873146 (growth habit), and SVU04_39821223 (maturity) are within 7.6 Mbp of each other. Other QTN hotspots occurred on Chromosome 8 (dry pod color and antioxidant content) and Chromosome 9 (dry pod color and maturity).

4. Discussion

Genotyping by Sequencing (GBS) showed high levels of SNP variations among cowpea samples, with 76K loci detected and almost 7K used in diversity analysis. The percentage of each SNP type in our study was 29.3, 29.0, 10.1, 12.2, 9.2, and 10.1% for [AG], [CT], [GT], [AT], [AC], and [CG], respectively. Transitions ([AG], [CT]) were more prevalent than transversions ([GT], [AT], [AC], and [CG]) in our cowpea panel.

Population structure analysis in our study differentiated the germplasm collection into two groups irrespective of geographical origin or taxa. This is in agreement with some previous results where diversity analysis of cowpea collections resulted in two or three major groups depending on whether breeding lines were included. For example, in their historical study, Herniter et al. [54] identified six subpopulations involved in the global spread of domesticated cowpea, using 368 accessions to represent worldwide diversity. Our study was similar in that one subpopulation deviated from the other five early during the domestication process, forming two groups at a higher level than among the six groups. Grain cowpeas and yardlong bean accessions grouped separately from each other in this study as in our previous one with both types [9].

Wild cowpea genotypes from V. unguiculata ssp. cylindrica and ssp. pubescence taxa clustered with both cowpea morphotypes, indicating a continued introgression amongst them and with ssp. unguiculata and ssp. sesquipedilis. However, the current study found that more variation existed within each subspecies rather than among subspecies groups (Table 1). This illustrates the availability of genetic resources to continue incorporation of favorable alleles from wild cowpea germplasm. Moreover, hybridization between wild and domesticated cowpea accessions can counteract narrowing genetic diversity precipitated by inbreeding a selection in vegetable or grain cowpea groups.

In terms of MTAs detected, the GWAS analysis revealed significant SNPs associated with all nine traits. Most of the significant SNPs were associated with the pod-length trait, for which there was high phenotypic variability, clearly indicating that multiple genetic mechanisms control this characteristic in cowpea. Despite the quantitative inheritance, we did note that one major genomic region was strongly associated with pod length. This region contained a cluster of significant SNPs and was located on Chromosome 1. A major QTL in the same region has been reported before following linkage mapping of bi-parental populations of yard-long bean [18,55]. However, our results also suggest there is another major region in Chromosome 10, and multiple minor regions distributed across the genome on Chromosomes 2, 3, 5, 6, 7, 8, 9, and 11. Single nucleotide polymorphisms on Chromosome 3 in a yardlong bean study [56] and Chromosome 8 in a cowpea domestication study [40] were also reported as associated with pod length, indicating potential minor QTL regions responsible for pod length in both groups.

Among the other traits measured, such as pod color, we found MTAs on Chromosomes 3, 5, 6, 8, and 9. These contained SNPs significantly associated with dry pod color. Purple pod color is dominant in cowpea and is controlled by at least two genes [57] and there is a strong correlation of pod tip color and seed coat color [32]. Our results indicate that pod color is indeed oligogenic and that significant SNPs are interspersed in genomic regions containing transcription factors associated with abiotic stress, such as zinc finger DNA-binding domains (Vigun03g106800.1) and AP2/B3-like transcriptional factors (Vigun03g107300.1) as well as senescence (Vigun03g405300.1) and plant development (Vigun03g406900.1).

Flower color is another qualitative important trait in cowpea, as it influences pollinator behavior [58] and is highly correlated to pod and sometimes seed coat color. All significant SNPs observed in this study for flower color were on SCONTIG_3 or Chromosomes 1, 3, 8, and 9. On Chromosome 8, a flavonoid biosynthesis gene (Vigun08g040200.1) was identified. This gene codes for Chalcone-flavanone isomerase family proteins that isomerize narigenin chalcone to naringenin, a primary precursor in flavonoid synthesis [59]. Flavanols are a group of flavonoids responsible for white or yellow floral colors upon UV absorption [60]. White and yellow are typical petal colors in cowpea, but purple is dominant over white [61]. Other genes associated with significant SNPs for flower color were related to plant signaling regulators or abiotic stress response. These included a heat shock transcription factor A2 (Vigun08g039300.1), thioredoxin superfamily proteins (Vigun01g178500.1, Vigun01g178600.1, and Vigun08g039200.1), and ovate proteins (Vigun01g178000.1 and Vigun01g178100.2) that regulate ovule development in various plant families [62,63].

Although seed color was not evaluated, the antioxidant content in cowpea seeds is a highly variable trait related to flavonoid pathway that also highly influences segmental colors found in the seed coat, especially around the hilum and ‘eye’ of the seed [31,64]. In our study, a total of 5 SNPs on Chromosomes 2, 4, 8, 9, and 10 were significantly associated with seed antioxidant content. Antioxidant levels associated with SNPs from Chromosomes 2 and 4 were located close to secondary metabolite catalytic enzyme genes including S-adenosyl-L-methionine-dependent methyltransferases (Vigun02g102800.1 and Vigun04g174900.1) and alpha-galactosidase 2 (Vigun04g174100.1).

Antioxidants are important components in the adaptive mechanism to oxidative stress in plants [65,66]. The SNP SVU08_36638094 on Chromosome 8 was identified in both the BLINK and FarmCPU models and was close to a genomic region populated with pectin lyase-like superfamily protein genes (Vigun08g202400.1, Vigun08g202800.1, Vigun08g202900.1, Vigun08g203000.1, Vigun08g203100.1, Vigun08g203200.1, and Vigun08g203300.1). Pectin lyases cleave glycosidic bonds in pectin molecules to produce pectin oligosaccharides, which have antioxidant properties [67]. Similarly, several defense related genes were also located near significant SNPs for antioxidant content such as leucine-rich repeat protein kinases (Vigun02g104600.1 and Vigun09g241400.1), aquaporin-like superfamily protein (Vigun04g174300.1), ethylene response factor 7 (Vigun08g202300.1), and allene oxide synthase (Vigun09g243000.1). Chromosome 8 also contained a major region of SNPs associated with seed protein content [68], which could be stacked in a breeding scheme to generate high-nutritional-value cowpea lines.

In terms of phenological traits and the evaluation of timing of developmental stages in cowpea, there were a few MTAs of interest. Namely, significant SNPs associated with pod maturity were found on Chromosomes 3, 4, 5, and 9. These were located in similar regions as significant MTAs for pod maturity detected by Andrade et al. [69] using the Cowpea iSelect Consortium SNP array, although exact SNP positions and genes differ. Key genes for pod maturity identified in the present study were related to a senescence, RING/U-box superfamily protein (Vigun04g067100.1), a cyclin family protein (Vigun05g135900.1); a plant development, late embryogenesis abundant protein (Vigun03g292500.1), an ARM repeat superfamily protein (Vigun04g066500.1), and a kinetochore protein (Vigun09g089700.1).

Growth habit in cowpea is predominantly erect, although some accessions are semi-prostrate or prostrate. Cowpea-producing regions across the globe have adopted a particular growth habit depending on end-user preferences and breeder choices [69]. Only one significant SNP was associated with cowpea growth habit in our study and was on Chromosome 4, but most genes around this SNP were of little interest as they were for uncharacterized or hypothetical proteins. The lack of more significant SNPs for growth habit could be due to the photoperiod of the testing site, the method of evaluating the trait, or reduced diversity of growth habit within the mini-core of 373 genotypes tested here. Growth habit QTL of high importance to yield and maturity date have been found in common bean, a New World relative of the cowpea, but one with more variability in determinacy and level of indeterminate vine growth [70,71].

One major yield component that was studied here, and which produced interesting MTAs, was seed weight. Significant SNPs associated with seed weight were distributed across Chromosome 3 (SVU03_46623413, SVU03_50505045, and SVU03_52946137) and Chromosome 4 (SVU04_16299944). The genes Vigun03g332100.1, Vigun03g332200.2, Vigun03g332400.1, and Vigun03g332500.1 all encode major facilitator superfamily proteins, broad spectrum substrate, and solute transporters in cells [72,73]. Other membrane transporter genes identified include the proton exporter H (+)-ATPase 11 (Vigun03g332000.4) and signal transduction genes calcium-dependent protein kinase 1 (Vigun03g331900.1), Ubiquitin carboxyl-terminal hydrolase-related protein (Vigun03g331800.1), and GRAS family transcription factor (Vigun03g285400.1). This may indicate that starch storage in seeds requires the action of cell metabolite transporters regulated by signal transduction genes and transcription factors.

Lucas et al. [74] identified 10 QTLs associated with 100-seed weight across nine mapping populations, most of which accounted for more than 30% of the phenotypic variance. Two of these QTLs were located on Chromosome 3 and Chromosome 4. Of these two QTLs, only the one located on Chromosome 4 was syntenic to soybean genome regions associated with seed size. The study by Lo et al. [33] identified 13 significant regions with candidate genes that included a cell wall protein, a phosphate transporter, a polycomb group protein, a histidine kinase 2, a WD repeat protein, and a delta (24)-sterol reductase.

Response to low phosphate and growth with rock phosphate had only one significant SNP each. The rock phosphate responsive gene was on Chromosome 11, located in a genomic region rich in abiotic stress response genes including AP2/B3-like transcriptional factor family protein (Vigun11g169600.1 and Vigun11g169700.1), P-loop, containing nucleoside triphosphate hydrolases superfamily protein (Vigun11g168200.1), Ypt/Rab-GAP domain of gyp1p superfamily protein (Vigun11g168300.1), and RWP-RK family transcription protein (Vigun11g168000.1). the SNP for low phosphate response was on Chromosome 10 and located along stress response genes such as Receptor-like protein kinase 1 (RPK1) (Vigun10g029500.2 and Vigun10g029600.1), and serine carboxypeptidase-like 20 (Vigun10g029300.1). The RWP-RK gene controls symbiotic root nodule development in soybean based on low nitrogen or Phytophthora sojae infection [75]. Another gene in soybean GmBBE-like43 coordinates root responses to both P deficiency and aluminum toxicity [76]. This may suggest a potential role in aluminum toxicity tolerance and response for the genes near the SNP associated with low phosphorus tolerance in this study.

5. Conclusions

Our work describes significant SNPs and candidate genes associated with nine very important traits in cowpea. All SNPs can be implemented in marker-assisted selection in any cowpea breeding program to accelerate the development of novel cowpea cultivars. Key genes contributing to variation of each trait have been identified and can be validated through functional genomics. The results underscore the importance of maintaining plant phenotyping databases at public institutions such as the USDA and strengthening collaborations to bolster data collection efforts, homogenizing data across sites, and maximizing the utility of limited phenotyping and genebank curation resources.

Finally, some co-localization of significant QTNs was observed in our study that could indicate genomic regions or genes regulating multiple traits pleiotropically or linked genes for multiple traits in tight LD. These were on Chromosome 3 for seed weight, flower color, maturity, and pod length; on Chromosome 4 for seed weight, antioxidant content, growth habit, and maturity; on Chromosome 8 for dry pod color and antioxidant content; and on Chromosome 9 for dry pod color and maturity.

The co-localization of SNPs associated with important domestication-related traits resemble the findings of Lo et al. [33], who observed co-localization and postulated some concerted molecular evolution. For example, SVU08_36154680 (dry pod color) and SVU08_36638094 (antioxidant content) MTAs are very close (0.4 Mbp) on Chromosome 8, while SVU09_16383249 (dry pod color) and SVU09_16383398 (maturity) were within 149 bp on Chromosome 9. Several previously identified SNPs associated with antioxidant content also co-localized with seed color [31]. Chromosomes 4 and 9 also harbored SNPs associated with cowpea flower color or flowering time [23,60,61] and might be relevant for comparisons. These genomic regions for antioxidant content or maturity date in cowpea could provide an opportunity for marker-assisted selection and multiple trait integration that can be prioritized for crop improvement.