Assessing the genetic diversity of cowpea [Vigna unguiculata (L.) Walp.] germplasm collections using phenotypic traits and SNP markers

Productivity of cowpea [Vigna unguiculata (L.) Walp] in sub-Sahara Africa is curtailed by a lack of farmer-preferred and improved cultivars and modern production technologies. The objectives of the study were to determine the extent of genetic diversity present among a collection of cowpea accessions from Zambia and Malawi using phenotypic traits and single nucleotide polymorphism (SNP) markers and, to select distinct and complementary parental lines for cultivar development. One hundred cowpea genotypes were evaluated for agronomic traits in two selected sites in Zambia, using a 10 × 10 alpha lattice design with two replications. Ninety-four of the test genotypes were profiled with 14,116 SNP markers. Number of pods plant− 1 (NPP), pod length (PDL), and number of seeds pod− 1 (NSP), were significantly (p < 0.05) affected by genotype × environment interaction effects. Genotypes such as CP411, CP421, CP645, CP732, Chimponongo, and MS1–8–1-4 exhibited higher grain yield of > 1200 kg/ha with excellent performance in yield components such as NSP, PDL, HSW and GYD. Grain yield had significant (p < 0.05) associations with NPP (r = 0.50), NSP (r = 0.46) and PDL (r = 0.42) useful for simultaneous selection for yield improvement in cowpea. The SNP markers revealed gene diversity and polymorphic information content of 0.22 and 0.17, respectively, showing that the tested cowpea accessions were genetically diverse. Test genotypes were classified into four genetic groups irrespective of source of collection allowing selection and subsequent crosses to develop breeding populations for cultivar development. Genotypes Bubebe, CP411, CP421, CP645, Chimponogo and MS1–8–1-4 were identified to be the most genetically divergent and high yielding making them ideal parental lines for breeding. This study provided a baseline information and identified promising cowpea genetic resources for effective breeding and systematic conservation.


Background
Cowpea [Vigna unguiculata (L.) Walp., 2n = 2x = 22] is a relatively low cost source of plant-derived protein, amino acids and essential nutrients globally. It is the main food staple supporting millions of people in sub-Sahara Africa (SSA) [1,2]. The grain protein content of cowpea is about 250 mg/g [3], which is comparable to that of soybeans [2]. In addition, cowpea grain contains essential nutrients such as iron (53.2 mg/kg), zinc (38.1 mg/kg), calcium (826 mg/kg) and magnesium (1915 mg/kg) [3]. Young and succulent leaves and pods of cowpea are used as cooked vegetable, while the grains are ground and processed into powder for making thick porridge, gravy or sometimes consumed as a boiled delicacy [4].
Cowpea is a key companion crop in mixed cropping systems useful to supressing weed infestation, enhancing soil fertility and reducing water evaporation [5]. Cowpea forms symbiosis with the root nodule bacterium, Rhizobium, and fixes 70 to 350 kg/ha of atmospheric nitrogen and some 40 to 80 kg of this is deposited into soils as a natural source of mineral nitrogen contributing to soil health [5]. Cowpea thrives under low soil fertility and dry-land growing conditions making it one of the most resilient legume crops suitable for the low input and water-limited production systems in SSA.
Global production of cowpea is estimated to be 6.5 million tons per annum on 14.5 million hectares of land [6]. The leading world producers of cowpea are Nigeria and Niger with five and three million hectares of production areas, respectively [7]. Cowpea is widely cultivated by small-scale farmers in southern African countries such as in Zambia, Zimbabwe, Malawi, Namibia, Mozambique and Botswana [8,9]. The mean grain yields of cowpea in SSA is between 100 to 599 kg/ ha which is far less than the potential yield of the crop reaching up to 3 t/ha elsewhere [8,10]. The yield gap is attributable to a lack of improved and high yielding cultivars, poor agronomic practices and an array of abiotic and biotic production constraints. Therefore, there is need to develop best performing, locally adapted and farmer-preferred cowpea varieties for sustainable production in the region.
The southern African countries including Namibia, Botswana, Zambia, Zimbabwe, Malawi, Mozambique and South Africa are believed to be the centres of diversity of cowpea where primitive and wild relatives are found [11]. Diverse cowpea germplasm collections are conserved in the Southern African Development Community (SADC) gene bank in Lusaka/Zambia. The country serves as Plant Genetic Resources Centre coordinating the works of some 16 National Plant Genetic Resources Centres (NPGRCs) in southern Africa [12]. Farmers in southern Africa widely grow unimproved landraces due to a lack of improved and locally adapted farmer-preferred cultivars. Landraces exhibit low yield potential, heterogeneous in flowering and maturity, poor processing quality, and low palatability and digestibility [13]. Low palatability and digestibility are adaptive traits against field and storage pests, traits resulted from repeated cycles of natural and artificial selection. The low palatability and digestibility of landraces reduce their utility for human consumption due to prolonged cooking time and reduced bioavailability of essential nutrients. Therefore, the cowpea genetic resources found in the region can be explored as a novel source of genetic variation for breeding programs.
A well-characterised crop genetic resource is a precondition for effective breeding and genetic conservation. Genetic diversity is assessed using phenotypic traits and molecular markers. Phenotypic characterisation in the target production environment enables identification and quantification of genetic variation for key qualitative and quantitative traits for ideotype breeding. Knowledge of phenotypic variation and traits relationship assist crop breeders to develop the most adaptive and productive cultivars [14]. The genetic diversity of cowpea for phenotypic traits is assessed using standard descriptors developed by the International Board for Plant Genetic Resource [15]. Key phenotypic traits include days to flowering, time to maturity, growth habit, flower colour, number of pod plant − 1 , pod length, number of seeds pod − 1 , seed colour, seed size, hundred seed weight and grain yield [15].
Various DNA markers such as the restriction fragment length polymorphism (RFLP), amplified fragment length polymorphism (AFLP), simple sequence repeat (SSR), random amplified polymorphic DNA (RAPD) and single nucleotide polymorphisms (SNP) have been used in cowpea genetic diversity analysis [16][17][18]. SNPs are markers of choice in genetic diversity analysis because they are widely distributed throughout the genome and their detection is amenable to automation [19]. In addition, SNP markers are increasingly time and cost efficient to genotype large populations with a relatively higher throughput [20]. SNP markers were applied in genetic diversity analysis of cowpea [16].
Cowpea is one of food security crops in Zambia widely cultivated in the eastern, southern and western regions. Hitherto, only seven cowpea varieties were released in the country that are relatively poor performers (< 700 kg/ha) and largely succumbed to emerging pests and diseases. The genetic diversity present among the germplasm collections conserved in the gene bank and landraces cultivated by smallholder farmers in Zambia can be explored for cowpea breeding and new cultivar deployment. Therefore, the objectives of the present study were to determine the extent of genetic diversity present among a collection of cowpea accessions from Zambia and Malawi using phenotypic traits and SNP markers, and to select distinct and complementary parental lines for cultivar development in Zambia.

Analysis of variance based quantitative phenotypic traits across locations
The combined analysis of variance revealed that the genotype × site interaction effects were significant (p < 0.05) for PDL, NPP and NSP (Table 1). DTF, DTM, PDL and NPP varied significantly (p < 0.05) between the two sites. The genotypes had varied flowering and maturity date as revealed by the significant (p < 0.05) genotypic effect. Similarly, there was significant (p < 0.05) genotype difference for PDL, NSP, HSW and GYD.

Variation based on qualitative phenotypic traits
There were significant differences (P < 0.00) among test genotypes for key qualitative traits (Additional file 1). For growth habit, 43 of the accessions were indeterminate, 39 determinate and 18 creeping types. Genotypes with predominantly upright growth type and short plant height were Bubebe, Namuseba, Msandile and MS1-8-1-4. Chimponongo and BBXSC13 had creeping growth type. Forty-nine accessions had brown and 21 black seed coat colour, while the rest of the genotypes had 12 purple-brown, 10 white and 8 red-brown. Based on leaf colour genotypes were assorted into light green (26 genotypes), light green (35) and dark green (39). Pod colour was variable varying from deep green (52 genotypes), light green (30) and purple (18). There were three classes of genotypes based on flower colour: 95 genotypes displayed violet flower, while four had yellow and one had white. Therefore, a combination of the assessed qualitative traits are useful markers for genotype selection in cowpea improvement programs.

Variance components and heritability of quantitative agronomic traits
Phenotypic coefficient of variation (PCV) values were higher than genotypic coefficient of variation (GCV) for all the traits ( Table 3). The GCV values ranged from 0 to 14.6%, while the PCV ranged from 0 to 21.56%. Larger discrepancies between GCV and PCV estimates were observed for all assessed traits. The genotypic variance accounted for ≥50% of the total variation for grain yield. Low heritability (≤ 30) estimates were recorded for days to maturity, hundred seed weight, number of seed pod − 1 Note: *, **,*** = Significance at 10, 5 and 1%, respectively; ns not significant, DF degrees of freedom, DTF days to flowering, DTM days to maturity, PDL pod length (cm), NPP number of pods per plant, NSP number of seeds per pod, SDS seed size, HSW hundred seed weight (g) and GYD grain yield in kg per hectare Note: CV coefficient of variation, LSD least significant difference, SE standard error, DTF days to flowering, DTM days to maturity, PDL pod length (UNIT?), NPP number of pods per plant, NSP number of seeds per pod, SDS seed size (mm), HSW hundred seed weight (g/100 seed), GYD grain yield in kg ha − 1 and pod length and number of pod plant − 1 . The heritability estimates for days to flowering and seed size were moderate (30-60%), while grain yield recorded heritability estimates above 60% that will enhance the response to selection and breeding gains. Genetic advance ranged from 0 to 20.58%. Seed size and days to flowering had moderate GA% (10-20%).

Principal component (PC) and bi-plot analyses
The first three PCs with Eigen-values greater than 1 accounted for 71.25% of the total variation exhibited by the assessed quantitative traits ( Table 5). The first principal component (

Genetic diversity and population structure
The tested SNP markers were moderately polymorphic with a mean PIC value of 0.17 for the whole population and 0.21 for the mutant lines ( Table 6). The mutant lines also exhibited the highest gene diversity (GD) with a mean of 0.26. The highest mean minor allele frequency was 0.22 observed among landraces. The whole population had high heterozygosity (0.30). Among the biological types, mutant lines exhibited the highest heterozygosity estimate with a mean of 0. 35 The structure analysis based on the Evanno method allocated the test genotypes into four main clusters with the highest value of ΔK that occurred at K = 4 (Fig. 2 a). Genotypes that scored < 0.80 were considered as pure line populations, while those that were < 0.80 as admixtures (Fig. 2 a). The model-based clustering using the 90 accessions showed the four admixture sub-populations (Fig. 2 c). Sub-population I was composed of 16 accessions (17.7%) that were sourced from Malawi and the University of Zambia. About 22 accessions (24.4%) were allocated in sub-population II and these genotypes were mainly acquired from Malawi, the National Gene Bank of Zambia and the University of Zambia. Sub-population Table 4 Correlation coefficients of grain yield and yield components among 100 cowpea genotypes evaluated at two locations in Zambia   (Fig. 2b). The analysis of molecular variance (AMOVA) showed a significant variation within populations. The within population variation accounted for 92% of the total variance ( Table 7). The variation detected among the populations was not significant and explained only 8% of the variability in the germplasm. The lack of genetic variation among the populations was confirmed by the low pair-wise genetic differentiation (Fst) values ranging between − 0.004 and 0.012 and inbreeding coefficient (Fis) of − 0.351 to − 0.362 (Table 8).

Combined analysis of phenotypic and genotypic data
The dendrograms for phenotypic and genotypic data each revealed three heterogeneous clusters among the genotypes. The phenotypic dendrograms showed that the first cluster comprised of genotypes from all sources,   while the second cluster comprised of elite lines from Malawi, University of Zambia and the National Gene Bank and the last cluster contained accessions from the National Gene Bank, Malawi and the University of Zambia (Fig. 3). The largest cluster in the genotypic dendrogram contained genotypes from the National Gene Bank, while the second largest consisted of mostly genotypes from Malawi (Fig. 4). The joint matrix revealed three similarly sized clusters among the genotypes (Fig. 5). The phenotypic and genotypic dendrograms were compared using the tanglegram and only a few genotypes (about 10%) maintained their positions (Fig. 6). The tanglegram comparison highlighted that both the positions and groupings of the genotypes were not consistent across the phenotypic and genotypic dendrograms. A relatively few genotypes maintained their position, i.e., about 30% maintained their groupings. Furthermore, the phenotype and genotype dissimilarity matrices exhibited a very low correlation (r = − 0.03) when subjected to the Mantel test (Additional files 2 and 4). In contrast, the genotype and phenotype dissimilarity matrices each exhibited strong correlations of r = 0.12, and r = 0.99 with the joint matrix, respectively (Additional files 5 and 6).

Genotypic variation and performance of test genotypes for key qualitative and quantitative traits
The present study evaluated the genetic diversity present among 100 diverse genotypes of cowpea germplasm collections from southern Africa using qualitative and quantitative phenotypic traits in two locations in Zambia. Further, high density SNP markers were used as a preliminary step to identify suitable and complementary parental lines for breeding. There were significant genotype × location interaction (Table 9) effect signifying that the tested germplasm were genetically diverse for selection and cultivar development targeting the test locations. Also, the interaction effect shows that the genotypes responded differently in the test environments which can facilitate identification of cowpea lines with specific or broad adaptation. Specific and broad adaptation have been identified and exploited in the Brazilian cowpea breeding programs based on genotype × location interaction analysis [21]. The interaction effect suggests that the test environments influence genotypic performance, which may confound genotype selection efforts by reducing the correlation between genotype and phenotypic expression [22].  In the present study, the assessed quantitative traits were affected by genotype × location interaction effect. Hence, there is intrinsic genetic variation influenced by the test locations necessitating multi environment evaluation for selection. Differential genotype response to environmental conditions during germplasm evaluation is attributable to the differences in genetic constitution among test genotypes and micro-environmental conditions [23]. In this study, the SCCI (mean yield of 832 kg ha − 1 ) site is high yielding environment compared with the GART (mean yield of 764 kg ha − 1 ) site probably due to the prevailing favourable environmental conditions such as better soil fertility and higher moisture levels in the former. Genotype phenology and biomass production exhibit environmental plasticity due to variable soil and climatic factors [24,25]. In the present study, some genotypes were high grain yielders (e.g. CP411 with 2197 kg ha − 1 ) and others were low yielders (e.g. ZM6680 with 87 kg ha − 1 . Quantitative traits are under the influence of polygenes. Hence, it is pertinent for genotype selections in multiple test environments to minimise environmental variance and to enhance selection gains [26,27]. Genotypes such as MS1-8-1-4, Msandile, BBXSC13, CP411, CP421, CP654, CP3413 and Bubebe that exhibited early to medium maturity are ideal candidates for drought tolerance breeding to offset the incessant droughts experienced in southern Africa. Early maturity is associated with drought escape [28] [8]. reported that farmers in southern Africa prefer cultivars with a short flowering period and maturity, valuable traits to evade the "hunger period". Highest number of seeds per pod (e.g. expressed by genotypes CP421 and Bubebe) is one of the factor affecting genotype responses based on their efficiency in growth resource utilisation and allocation. This could also be contributed to increased length of the pods by test genotypes [29]. Seed weight is directly associated with seed size and it is recommended to be used as an indirect selection criterion to maximise grain yield response in cowpea [10] [8]. reported a high yield potential of cowpea genotypes that can reach up to 3 t/ha. The yield level recorded in the present study by the landraces was relatively less. This   could be the low yield potential of landraces grown by most farmers in SSA. In the region, landraces are continuously cultivated because they possess farmerpreferred quality traits and their ability to adapt under variable stress conditions due to their genetic diversity and plasticity [30,31].
In the present study qualitative traits such as seed coat, pod and leaf colour were more important traits for selection. These traits affect the market value of cowpea in Africa given that farmer and consumer preference are based on these attributes. Seed coat colour is often associated with processing quality (e.g. cooking time) and farmers deliberately select varieties that have shorter cooking time [32]. The inheritance of seed coat colour is governed by few major genes that will enhance selection progress during cultivar development [33]. In this study, the genotypes Bubebe and Msandile, with predominantly light-green leaves exhibited determinate growth habit in comparison with BBXSC13 and Chimponongo that had dark green leaves and creeping growth habit [34]. reported that cowpea cultivars with a determinate growth type were more drought tolerant compared to the indeterminate types [35]. reported that indeterminate varieties of cowpea attained higher productivity due to their prolonged maturity and photosynthesis efficiency. Therefore, in order to promote sustainable production and productivity and enhanced adoption of improved cowpea cultivars, breeding programs should incorporate farmer-and market-preferred attributes in the newly developed cultivars.

Variance, heritability and genetic advance
In this study, the heritability estimate for grain yield was high (64%), suggesting that the grain yield achieved by the accessions was highly repeatable ensuring genetic improvement through selection. The high heritability value for grain yield corroborates with the findings of [2] but lower than a heritability value of 97% reported by [29]. Genetic advance is directly related with yield gains achievable via selection. High estimates of genetic advance (e.g. for HSW and SDS) and high heritability indicate that selection would result in foreseeable genetic improvement [36,37]. The large discrepancies values of PCV compared to GCV in this study, suggests that trait expression was also influenced by environment factors in addition to genetic effects, which was also confirmed by the significant location main effects in the ANOVA ( Table 2).

Associations of quantitative traits
The relationships among yield and yield components are critical in devising a selection strategy. Selection of one trait may amplify or negatively affect performance in the other traits. The high contribution and strong association of PDL, NSP, GYD, and NPP to PC1 as well as DTF and DTM with the PC2 indicated that these traits were highly discriminatory explaining the variation among the genotypes [38,39]. found that traits such as NPP and GYD in cowpea were associated with PC1 showing the importance of agronomic traits in cowpea evaluation corroborating with the findings of the present study. The genotype-trait biplot enables visual and simultaneous selection of genotypes for multiple traits. There was strong correlations between PDL, NPP, NSP and grain yield indicating their positive impact on genotype performance. Previous reports identified these traits being important yield-influencing attributes [29,39]. Entries such as E10 (CP 411), E71 (LT16-7-2-5), E13 (CP421) and E20 (CP645) scored greater grain yield response and yield-influencing traits suggesting their utility in variety improvement for yield gains and breeding population development. Entries such as E10 (CP411), E20 (CP645), E13 (CP421) and E58 (Sundan1) are selected with desirable NSP, GYD, PDL and DTF, respectively.

Population structure and genetic parameters
Genetic analysis using SNP markers delineated the test populations in to four genetic groups. This demarcation was irrespective of source of collection, suggesting that geographical sources of collection are not the sole factor for classification of cowpea genotypes. The germplasm exhibited high heterozygosity of 0.30, which showed that alternative alleles were represented in this population. This could be attributed to genetic exchange between and among cowpea breeding programs in sub-Sahara Africa, particularly in southern Africa. Heterozygosity estimates were highest for mutant lines, followed by elite lines and landraces. This suggests that mutant lines exhibited more variability possibly due to continuous gene segregation and other chromosomal aberrations. Landrace varieties are not subjected to artificial selection by breeders. Hence, landraces of self-pollinating crops including cowpea consists of heterogeneous populations due to genetic variation until complete homozygosity is reached through selection [40]. Mutant and elite lines are the result of several cycles of deliberate selection and inbreeding with expectedly high level of homozygosity at multiple loci. The present mutant lines were grouped together with the mother stocks or parental lines. This shows that the mutation events were minor and the genetic background largely remained the same for most of the assessed traits in the mutants lines compared with the elite parental lines. Malawi and Zambia are in close geographical proximity hence germplasm exchange between the two countries cannot be ruled out. Trait preference for the farmers and the market in the region may not be significantly different leading to the overlap of cowpea genetic resources in these agro-ecologies. This has partly disallowed the population structure analysis without distinguishing the genotypes based on geographical sources agreeing to the report of [18]. Exchange of genetic resources is key for plant-breeding research and cultivar development, which are dependent on wider genetic bases [10]. The PIC and GD values were essential to assess genetic diversity within the whole population, and subpopulations for identification selection of divergent parental lines useful for breeding. The moderate PIC and GD values exhibited by the whole population could be attributed to the inherent nature of self-pollination in cowpea, which limits genetic diversity. Self-pollinating species often lack genetic diversity due to higher levels Hundred seed weight HSW Weight of one hundred randomly selected seeds of a genotype measured in grams 8 Grain yield GYD The average grain yield per plot and converted into kg ha −1 using the formula given above.
of homozygosity at multiple loci. The within-population diversity in self-pollinating species is relatively low but varies among populations [41,42]. In contrast, mutants lines exhibited higher Ho, PIC and GD values compared to the elite lines or landraces. Naturally, landraces of self-pollinating species would be expected to exhibit higher levels of homozygosity compared to mutant and elite lines that were subjected to deliberate mutagenesis and crossing during breeding. On the hand, mutant lines exhibited higher heterozygosity and PIC values than elite lines, which could be due induced random mutations. Mutation breeding is a well-established technique to widen genetic diversity [43]. However, the mutant and elite lines displayed difference in growth pattern. This may be attributed to genetic differences in plant architecture in terms of growth habit and maturity period, among others. For example, the elite line Bubebe and the mutant line MS1-8-1-4, have short maturity period and determinate growth type. Hence, inclusion of landraces such as Chimponongo and Kapita with late maturity and creeping growth habit would be recommended to increase genetic variation and to enhance genetic gain through selection. This is consistent with the findings of other cowpea researches who indicated that architecture of the crop results in genetic diversification [9]. The mean Fst values recorded in the present study showed low genetic differentiation among the test populations (max Fst = − 0.012), which could be attributed to the self-pollination nature of cowpea. This shows that the sub-groups account for a small proportion of the total variance exhibited by the whole population and selection would be more efficient between genotypes rather than between sub-groups. However, the negative Fis values indicate that individuals accounted for a large proportion of the total variance due to the presence of heterozygotes that could be useful in future selection programs. This could due to high heterozygosity values estimated for the whole germplasm underpinned by mutants and elite lines that have undergone artificial selection, and landraces that may be exhibiting segregation. Individual selection will enhance genetic advance compared to selection within a sub-group. The phenotypic and genotypic distance matrices exhibited low and non-significant correlation, which confirm that phenotypic and genotypic matrices were independent but complementary. The discordance between phenotypic and genotypic matrices is partly attributed to the significant environmental effect on the variable expression of the assessed phenotypic traits. Phenotypic performance is the result of the genotype, environment and genotype-by-environment interaction effects [44]. In addition, the inconsistences between genotypic and phenotypic matrices could be due to the fact that genotyping may be able to capture subtle genetic mutations across the whole genome that may not be expressed phenotypically. Other studies reported inconsistences between phenotypic and genotypic matrices in common beans (Phaseolus spp) [45]), yam (Dioscorea spp) [46] and Brassica spp. [46]. Due to the expected discord between genotypic and phenotypic matrices, the use of a joint matrix derived from combined phenotypic and genotypic matrices has been recommended to increase precision [47,48]. The strong correlations exhibited by phenotypic and genotypic matrices each with the joint matrix show that each were derivative of the joint matrix and can be used for increased precision as they do not overlap. An increase of more than 150% in precision was reported while calculating dissimilarity distances using a joint matrix compared to the phenotypic matrix [49].

Conclusions
Phenotypic analysis using qualitative and quantitative traits and genotyping using high density SNP markers revealed the presence of significant variation among 100 cowpea germplasm collections of southern Africa. Trait association analysis revealed significant correlation between NPP, NSP, PDL and GYD that could allow direct selection to improve grain yield. The SNP markers used in the study were able to deduce genetic variation among the tested cowpea populations. The largest proportion of variation was attributable to individual genotype differences that is essential for improving grain yield by crossing lines from different divergent populations. Test genotypes were classified in to four genetic groups irrespective of source of collection allowing selection for subsequent cross combinations to develop breeding populations for cultivar development. Genotypes Bubebe, CP411, CP421, CP645, Chimponogo and MS1-8-1-4 were identified being the most genetically divergent and high yielding making them ideal parental lines for breeding. This study provided a baseline genetic profile and identified promising cowpea genetic resources for effective breeding and systematic conservation.

Plant materials
The study evaluated 90 cowpea genotypes acquired from different sources (Supplementary Table 1). The germplasm included 61 elites lines from Malawi, International Institute for Tropical Agriculture (IITA), Nigeria, The National Gene Bank of Zambia and the University of Zambia, 15 landraces collected from smallholder farmers in Zambia and 14 mutant lines, which were derived from three parental lines; Lutembwe, Bubebe and Msandile. The accessions from IITA and the released cultivars were used as standard checks.

Description of the study sites
The 100 genotypes were field evaluated during the 2017/ 2018 main crop season at the following two sites: the Seed Control and Certification Institute (SCCI) in Chilanga and Golden Valley Agricultural Research Trust (GART) in Chisamba/Zambia. The SCCI site is situated at a latitude of 15 o 32′S and a longitude of 28 o 11'E with an altitude of 1206 m above sea level. The total mean annual rainfall at the SSCI site is 1092 mm, while the mean daily minimum and maximum temperatures were 12°C and 26°C, respectively. The GART site is situated at a latitude of 14 o 96′S and a longitude of 28 o 10'E and an altitude of 1103 m above sea level. The GART site receives a total mean annual rainfall of 884 mm with mean daily minimum and maximum temperatures of 10°C and 30°C, respectively. The soils at both sites are classified as Haplustalf clays with pH of 5.8 and 5.2 at SCCI and GART, respectively [50].

Trial design, field planting and management, and data collection
The experiments were laid out in a 10 × 10 alpha lattice design with two replications. Each genotype was sown in a plot with two rows of 5 m long. The plot area was 3.75m 2 . The inter-row and intra-row spacings were 75 and 45 cm, respectively. Two seeds were sown per station at a depth of 2 cm and later thinned to one plant two weeks after emergence. Basal fertiliser (N: P: K), containing 20% nitrogen, 10% phosphorus and 20% potassium, was applied at a rate of 200 kg ha − 1 prior to planting. All other agronomic practices for cowpea production were followed as recommended for Zambia [50]. The crops were grown under rain-fed conditions and both sites received an annual rainfall of 850 mm during the study.

Data collection
Data was collected from six qualitative and eight quantitative traits following the descriptors of the [15,51]. The list of traits and details of data collection and units are provided in Table 2. Grain yield was determined in kg ha − 1 based on the following formula: plot weight plot area x 100 − 14 100 − mc Â 10; 000 where; mc is moisture content measured at harvesting, 14% is standard constant moisture content for legumes [51] and 10,000 is a conversion factor for a hectare.

Data analysis
The frequency of test genotypes displaying the assessed qualitative traits were summarised and statistical significant tests conducted using the cross tabulation procedure with the Statistical Package for the Social Sciences (SPSS) version 24 [52]. The quantitative data was subjected to analysis of variance (ANOVA) using the alphalattice procedure in GenStat® version 18 [53]. A combined analysis of variance was conducted after detecting significant differences among tested genotypes in each location. The following linear model was used for the combined analysis of variance: βijk = μ + Gi + Ej + Gi * Ej + Ei (rk)(b) + εijk, where; βijk = observed response; μ = grand mean Gi = the effect of i th genotype; Ej = the effect of j th location, Gi * Ej = the genotype x location interaction effect; Ej(rk)(b) = error associated with k th replication in blocks in the j th location and εijk = experimental error. The blocks within replications were considered as random factor, while genotypes and locations were fixed factors. Trait means of test genotypes were separated using the Fischers Unprotected LSD at p ≤ 0.05 significance level. Genotypic, genotype by location interaction and phenotypic variances were computed from the excepted mean squares of the analysis of variance as follows; σ 2 g ¼ ðmsg − mseÞ lr ; σ 2 gl ¼ ðmsgl − mseÞ r ; σ 2 p ¼ σ 2 g þ σ 2 e þσ 2 gl, where; σ 2 g = genotypic variance, σ 2 gl = genotype by location interaction variance, σ 2 : p = phenotypic variance, msg = mean square of genotype, mse = mean square of error, l = number of location and r = number of replication. Variances below zero were adjusted to zero according to [54]. Heritability in broad sense (H 2 ) was computed according to [55]; H 2 ¼ σ 2 g σ 2 p x 100 where; σ 2 g is genotypic variance and σ 2 p is phenotypic variance. Heritability was categorized as low (0-0.30), moderate (0.30-0.60) and high (> 0.60) following [56]. A covariance analysis was performed to calculate coefficient of variations. The genotypic coefficient of variation (GCV) and phenotypic coefficient of variation (PCV) expressed in percent were computed as described by [57] as follows: GCV ¼ ð ffiffiffi ffi σ 2 p g χ Þ Â 100; PCV ¼ ð ffiffiffi ffi σ 2 p p χ Þ Â 100 , were σ 2 g = genotypic variance, σ 2 p = phenotypic variance, χ = grand mean. Genetic advance was calculated following [58] as follows: GA = (k) (σp) (h 2 ), where, GA = Genetic advance; k = selection differential at 5% selection intensity; σp = phenotypic standard deviation; h 2 = broad sense heritability; Genetic advance as a per cent of mean (GAM) was computed following [59]: GAM ¼ ½ GA χ Â 100, where, GA = (%) = Genetic advance as a per cent of mean; GA = Genetic advance; χ = Grand mean. Genetic advance as a per cent of mean was classified and rated based on the scales given by [60] as low (< 10%), moderate (10-20%) and high (> 20%).
The magnitude of traits relationship was determined using Pearson's correlation coefficients (r) using the SPSS version 24 [52]. Principal component analysis (PCA) was performed using the same software to examine the number principal components and trait associations. The principal components (PCs) with Eigen-values ≥1.0 were considered to explain the variation in phenotypic traits among the genotypes. PCA biplots were constructed in GenStat [53] to depict the relationships among the studied genotypes and traits.

DNA isolation and genotyping
Ten seeds of each cowpea genotype were planted in a plastic pot. The seedlings were allowed to grow to the three-leaf stage before fresh leaves were harvested for DNA extraction. Leaves were sampled from each genotype for DNA extraction. Fifty milligrams of total genomic DNA was extracted from the well-developed trifoliate leaves with the NucleoSpin plant II kit (Macherrey-Nagel, Duren, Germany) using the Lysis Buffer 1 (based on the CTAB method) according to the manufacturer's procedures. The DNA concentration of each sample was measured using a NanoDrop 1000 (Invitrogen, California, USA). For verifying DNA integrity, 2 μL of DNA were subjected to gel electrophoresis on 1.0% (w/v) agarose gel, stained with ethidium bromide. Subsequently, 40 μL of a 50 ng/μL DNA of each sample were genotyped with Illumina Cowpea iSelect Consortium Array using Diversity Arrays Technology (DArT) markers. In total, 94 cowpea genotypes were genotyped by the genotyping by sequencing (GBS) technology as described by [61] with 20,000 DArT markers. The markers were integrated into a linkage map by inferring marker order position from the consensus Dart map. Genotyping of the materials was carried out at the Biosciences eastern and central Africa-International Livestock Research Institute (BecA-ILRI) in Kenya.

Data analyses SNP filtering
For quality control, DArTseq SNP derived markers were filtered to remove bad SNPs and genotypes using the software's PLINK 1.9 in MS window and R statistical package. Markers and genotypes with > 20% missing data were eliminated. Rare SNPs with < 5% minor allele frequencies were also pruned from the data. Only 14,116 informative SNP markers and 90 genotypes were used for analysis after removal of 4240 SNPs and four genotypes. The four genotypes, CP1, CP2, CP479, and CP2223 were removed due to extreme heterozygosity (< 90%), duplication or high levels of missing data (> 20%). The chromosomal coverage of the 14,116 SNPs was presented in Supplementary Fig. 3.