Genetic Architecture of Salt Tolerance in Cowpea (Vigna unguiculata (L.) Walp.) at Seedling Stage Using a Whole Genome Resequencing Approach

Cowpea (Vigna unguiculata (L.) Walp.) is a diploid legume crop used for human consumption, feed for livestock, and cover crops. Earlier reports have shown that salinity has been a growing threat to cowpea cultivation. The objectives of this study were to conduct a genome-wide association study (GWAS) to identify SNP markers and to investigate candidate genes for salt tolerance in cowpea. A total of 331 cowpea genotypes were evaluated for salt tolerance by supplying a solution of 200 mM NaCl in our previous work. The cowpea panel was genotyped using a whole genome resequencing approach, generating 14,465,516 SNPs. Moreover, 5,884,299 SNPs were used after SNP filtering. GWAS was conducted on a total of 296 cowpea genotypes that have high-quality SNPs. BLINK was used for conducting GWAS. Results showed (1) a strong GWAS peak on an 890-bk region of chromosome 2 for leaf SPAD chlorophyll under salt stress in cowpea and harboring a significant cluster of nicotinamide adenine dinucleotide (NAD) dependent epimerase/dehydratase genes such as Vigun02g128900.1, Vigun02g129000.1, Vigun02g129100.1, Vigun02g129200.1, and Vigun02g129500.1; (2) two GWAS peaks associated with relative tolerance index for chlorophyll were identified on chromosomes 1 and 2. The peak on chromosome 1 was defined by a cluster of 10 significant SNPs mapped on a 5 kb region and was located in the vicinity of Vigun01g086000.1, encoding for a GATA transcription factor. The GWAS peak on chromosome 2 was defined by a cluster of 53 significant SNPs and mapped on a 68 bk region of chromosome 2, and (3) the highest GWAS peak was identified on chromosome 3, and this locus was associated with leaf score injury. This peak was within the structure of a potassium channel gene (Vigun03g144700.1). To the best of our knowledge, this is one the earliest reports on the salt tolerance study of cowpea using whole genome resequencing data.


Introduction
Salt stress is a growing threat affecting cowpea production worldwide [1,2].Previous studies have shown salt stress resulting in significant crop yield losses [3][4][5].Soil salinity affects over 830 million hectares of croplands worldwide [6][7][8].In the U.S., soil salinity has been reported on over 19.6 million hectares of croplands.The cost associated with soil salinity is estimated to be 12 billion USD annually [9][10][11].In addition, irrigation water from rivers in semi-arid croplands can rapidly increase soil degradation through a continuous accumulation of salt [12].
Effects of salinity have been shown to be increasingly more insidious in semi-arid areas where cowpea cultivation is prevalent [13].For example, cowpea grown in the Coachella Valley of California has been negatively affected by salinity [14].Salinity can lead to significant plant phytological and morphological damage to cowpea plants [15].At the seedling stage, salinity reduces the photosynthetic activity of cowpea plants, where excess of Na + in plant tissue can lead to plant death [16,17].The negative impacts of salinity on crops are more significant in semi-arid and arid regions where cowpea cultivation is prevalent, which could refrain growers from expanding areas for cowpea production [18].Salt stress reduced cowpea plant vigor and yield-reducing for cowpea [19].Salt stress can also cause a significant decrease in yield for cowpea grown on calcareous soils [20].The seedling stage has been found to be the most vulnerable stage to salt stress in cowpea [8].Therefore, evaluating salt tolerance at the seedling stage and mapping QTL alleles for salt tolerance at the seedling stage are needed.
DNA sequencing by the next-generation sequencing technologies (van Dijk et al., 2014), using improved double digest Restriction Site Associated DNA sequencing (ddRADseq) [21] and Genotyping by Sequencing (GBS) [22] is used for single nucleotide polymorphism (SNP) discovery and genotyping.GBS and ddRADseq are affordable ways to develop SNPs and SNP genotyping, establish genetic maps, and map QTLs [22].However, both GBS and ddRADseq have disadvantages, such as untargeted and large amounts of missing data across given taxonomic samples due to their reduced genome sequencing.With the decrease in DNA sequencing cost, whole-genome sequencing (WGS) and whole-genome resequencing (WGR) provide opportunities to develop a large number of SNPs [23].WGS can generate whole-genome assemblies (genome sequences) for any plant with and without genome information through de novo assembly [22].WGR provides a high-resolution, base-by-base view of the genome.In addition, the WGR is widely used to detect single nucleotide variants (SNVs), including SNPs, insertions and deletions (InDels), structural variants (S.V.s), and copy number variation (CNV).It also allows the examination of SNVs, InDels, S.V.s, and CNVs in the genome's coding and non-coding regions with reliable sequence coverage and coverage uniformity at whole-genome level (http://www.illumina.com/techniques/sequencing/dna-sequencing/whole-genome-sequencing.html(accessed on 3 March 2020)).Therefore, we will use the WGR technology for genotyping and SNP discovery for cowpea salt tolerance in this study.
We have done several experiments for salt tolerance in cowpea at the germination and seedling stages [24].In one of the previous experiments, 155 cowpea genotypes as an association panel were phenotyped for foliar injury, plant height, and fresh and dry shoot weight under 0 mM and 200 mM NaCl (without salt stress and salt treatment) conditions and genotyped using 1049 SNPs postulated from genotype by sequencing (GBS).Association study in this panel showed three SNP markers, Scaffold87490_622, Scaffold87490_630, and C35017374_128, associated with salt tolerance at the germination stage, and seven SNP markers, Scaffold93827_270, Scaffold68489_600, Scaffold87490_633, Scaffold87490_640, Scaf-fold82042_3387, C35069468_1916, and Scaffold93942_1089, associated with salt tolerance at seedling stage [24].
In this study, we conduct GWAS with 331 cowpea genotypes as an association panel for salt tolerance at the early seedling stage to identify SNP markers associated with salt tolerance in cowpea using WGR technology.The phenotypic data of salt tolerance in the 331 cowpea genotypes have been published [25].Here, we will report SNP markers and candidate genes for salt tolerance in cowpea.

Leaf SPAD Chlorophyll under Salt Stress
Results indicated that a total of 65 SNPs were significantly associated with leaf SPAD chlorophyll under salt stress in cowpea (Figures 1 and 2).These SNPs were located on chromosomes 1 and 2. Chromosome 1 harbored a total of 9 significant SNPs, whereas chromosome 2 had a total of 56 significant SNPs (Table S1).LOD (−log 10 (p-value)) values varied from 7.53 to 10.68.The first locus that was identified to be associated with leaf SPAD chlorophyll under salt stress was defined by a cluster of significant SNPs mapped on a 3 kb region of chromosome 1.The second locus that was found to be associated with leaf SPAD chlorophyll under salt stress was defined by a group of significant SNPs mapped on an 890 kb region of chromosome 2.The significant SNPs that were found on chromosome 1 were Vu01_24245081 (LOD = 7.57), Vu01_24246312 (LOD = 8.00), Vu01_24246319 (LOD = 8.00), Vu01_24246550 (LOD = 7.76), Vu01_24246587 (LOD = 7.94), Vu01_24246822 (LOD = 8.27), Vu01_24246905 (LOD = 8.08), Vu01_24246981 (LOD = 8.07), and Vu01_24248242 (LOD = 7.85) (Figure 1).The SNP that was closest to an annotated gene, Vigun01g086000.1, was Vu01_24245081.Vigun01g086000.1 encodes for the GATA transcription factor whose predicted tertiary structure is shown in Figure 1.
Int. J. Mol.Sci.2023, 24, x FOR PEER REVIEW 3 of 21 331 cowpea genotypes have been published [25].Here, we will report SNP markers and candidate genes for salt tolerance in cowpea.

Leaf SPAD Chlorophyll under Salt Stress
Results indicated that a total of 65 SNPs were significantly associated with leaf SPAD chlorophyll under salt stress in cowpea (Figures 1 and 2).These SNPs were located on chromosomes 1 and 2. Chromosome 1 harbored a total of 9 significant SNPs, whereas chromosome 2 had a total of 56 significant SNPs (Table S1).LOD (−log10(p-value)) values varied from 7.53 to 10.68.The first locus that was identified to be associated with leaf SPAD chlorophyll under salt stress was defined by a cluster of significant SNPs mapped on a 3 kb region of chromosome 1.The second locus that was found to be associated with leaf SPAD chlorophyll under salt stress was defined by a group of significant SNPs mapped on an 890 kb region of chromosome 2.The significant SNPs that were found on chromosome 1 were Vu01_24245081 (LOD = 7.57), Vu01_24246312 (LOD = 8.00), Vu01_24246319 (LOD = 8.00), Vu01_24246550 (LOD = 7.76), Vu01_24246587 (LOD = 7.94), Vu01_24246822 (LOD = 8.27), Vu01_24246905 (LOD = 8.08), Vu01_24246981 (LOD = 8.07), and Vu01_24248242 (LOD = 7.85) (Figure 1).The SNP that was closest to an annotated gene, Vigun01g086000.1, was Vu01_24245081.Vigun01g086000.1 encodes for the GATA transcription factor whose predicted tertiary structure is shown in Figure 1.The second locus, defined by an 890 kb region of chromosome 2 harbored nine annotated genes.Table S1 shows that the SNPs with the highest LOD (−log 10 (p-value)) were
black and grey dots represent the SNPs.The x-axis is the chromosome number, and the y-axis is the LOD or −log10 of the p-value.The horizontal red and blue bars are two different LOD thresholds.Below the Manhattan plot are gene IDs from phytozome v.13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020)) corresponding to the significant locus on chromosome 2. Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast.cgi(accessed on 3 March 2020)).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 3 March 2020)) and presented on the left-hand side in the above figure.

Protein Homologs and Gene Ontology
Protein homolog search was investigated for the candidate genes with functional annotations that could be linked to tolerance to plant abiotic stress.In this study, the search was carried out across the genomes of legumes such as soybean, common bean, and Medicago.Proteins that have a similarity >90% with the query were taken into account.In order to estimate the number of copies of each candidate gene for cowpea, a search was conducted within the cowpea genome.For the candidate genes associated with leaf SPAD chlorophyll under salt stress, multiple copies of Vigun02g128900.1,Vigun02g129000.1, and Vigun02g129300.1 within the cowpea genome (Table 2).The candidate genes Vigun01g086000.1,Vigun02g128700.1,Vigun02g129100.1,Vigun02g129200.1,Vigun02g129400.1, and Vigun02g129500.1 had one to three copies within the cowpea genome.The number of protein homologs was highest within the soybean genome, whereas it was lowest within the Medicago genome (Table 2).For the candidate genes associated with relative tolerance index for chlorophyll, a large number of copies of Vigun10g104200.1,Vigun10g104300.1, and Vigun10g104400.1 were found within the cowpea genome.The candidate gene Vigun04g193500.1 was unique within the cowpea genome.The candidate genes Vigun01g086000.1,Vigun02g129000.1,Vigun02g129100.1,Vigun02g129200.1,Vi-gun02g129300.1,Vigun02g129400.1,Vigun03g118000.1, and Vigun10g093500.1 had one to four copies within the cowpea genome.Overall, the number of homologs between common bean and cowpea was very close.Among the four legume species compared in this study, the soybean genome had the largest number of copies.For the candidate genes associated with leaf injury score, the number of gene duplications is less significant compared to other traits.The candidate genes Vigun01g086000.1,Vigun03g144700.1,Vigun04g025900.1, and Vigun04g193700.1 were unique within the cowpea genome.The candidate genes Vi-gun04g193700.1,Vigun02g129000.1,Vigun02g129200.1,Vigun03g135800.1,Vigun03g136300.1,Vigun03g149400.1, and Vigun04g054000.1 had one to four copies within the cowpea genome.Vigun03g135800.1 seemed to be abundant within the soybean, common bean, and Medicago genomes.However, only one of the common bean genomes had a single copy of Vigun01g086000.1.

Overlapping SNPs and Functional Annotations
The number of overlapping SNPs between traits was visualized using a Venn diagram (Figure 11A).On the Venn diagram, the significant SNPs associated with leaf SPAD chlorophyll under salt stress, relative tolerance index for chlorophyll, and leaf score injury were represented using solid green, blue, and pink circles, respectively (Figure 11A).The number of SNPs associated with leaf SPAD chlorophyll under salt stress, relative tolerance index for chlorophyll, and leaf score injury were 65, 60, and 1667, respectively.
A total of 19 SNPs overlapped between leaf SPAD chlorophyll under salt stress, relative tolerance index for chlorophyll, and leaf score injury, as shown in Figure 11A, suggesting that there could be a common genetic mechanism controlling these traits.The number of common SNPs between leaf SPAD chlorophyll under salt stress and tolerance index for chlorophyll was 3. The number of overlapping SNPs between the relative tolerance index for chlorophyll and leaf injury score was 4. The number of shared SNPs between leaf SPAD chlorophyll under salt stress and leaf injury score was 30.These results provided strong evidence of the interdependency between these traits at the genetic level.
Overlapping functional annotations between candidate genes associated with leaf SPAD chlorophyll under salt stress, relative tolerance index for chlorophyll, and leaf injury score were also visualized using a Venn diagram (Figure 11B).Duplicated functional annotation names were removed, and only the number of unique names was displayed on the Venn diagram.Color coding was similar to Figure 11A.The three traits investigated for salt tolerance showed a common functional annotation, supporting the evidence of the potential common genetic mechanism controlling these traits (Figure 11).In addition, a common functional annotation was identified for the candidate genes associated with leaf SPAD chlorophyll under salt stress and relative tolerance index for salt stress.No common functional annotation was found between the candidate genes associated with leaf SPAD chlorophyll and leaf injury score under salt stress.Similar results were found in the candidate genes associated with relative tolerance for chlorophyll and leaf injury score under salt stress.A total of 19 SNPs overlapped between leaf SPAD chlorophyll under salt stress, relative tolerance index for chlorophyll, and leaf score injury, as shown in Figure 11A, suggesting that there could be a common genetic mechanism controlling these traits.The number of common SNPs between leaf SPAD chlorophyll under salt stress and tolerance index for chlorophyll was 3. The number of overlapping SNPs between the relative tolerance index for chlorophyll and leaf injury score was 4. The number of shared SNPs between leaf SPAD chlorophyll under salt stress and leaf injury score was 30.These results provided strong evidence of the interdependency between these traits at the genetic level.
Overlapping functional annotations between candidate genes associated with leaf SPAD chlorophyll under salt stress, relative tolerance index for chlorophyll, and leaf injury score were also visualized using a Venn diagram (Figure 11B).Duplicated functional annotation names were removed, and only the number of unique names was displayed on the Venn diagram.Color coding was similar to Figure 11A.The three traits investigated for salt tolerance showed a common functional annotation, supporting the evidence of the potential common genetic mechanism controlling these traits (Figure 11).In addition, a common functional annotation was identified for the candidate genes associated with leaf SPAD chlorophyll under salt stress and relative tolerance index for salt stress.No common functional annotation was found between the candidate genes associated with leaf SPAD chlorophyll and leaf injury score under salt stress.Similar results were found in the candidate genes associated with relative tolerance for chlorophyll and leaf injury score under salt stress.

Discussion
Whole genome resequencing has been more and more popular in plant genetic-related studies.It allows for the discovery of a large number of SNPs that can be used in GWAS.Thanks to the large number of SNPs, the likelihood of discovering good candidate genes is higher [26,27].This study was one of the earliest reports in cowpea using a whole genome resequencing data to conduct GWAS for salt tolerance in cowpea.Whole genome resequencing provided a total of 14,465,516 SNPs.GWAS was conducted using a total of 5,884,299 filtered and high-quality SNPs.
In this study, a total of 65, 60, and 1667 SNPs were found to be significantly associated with leaf SPAD chlorophyll under salt stress, relative tolerance index for chlorophyll, and leaf score injury, respectively.The first reported molecular markers associated with salt tolerance in cowpea were Scaffold87490_622, Scaffold87490_630, C35017374_128, Scaf-fold93827_270, Scaffold68489_600, Scaffold87490_633, Scaffold87490_640, Scaffold82042_ 3387, C35069468_1916, and Scaffold93942_1089 [24].These are SNP markers that were obtained from genotyping-by-sequencing.The sequence that contains these SNPs was realigned to the cowpea genome to find whether they are mapped in the vicinity of the SNPs reported in this paper.The SNPs Scaffold87490_622, Scaffold87490_633, and C35069468_1916 were found at about 10 kb downstream of Vu01_24246905, indicating that the results from this study are complementary with the GBS study that was previously reported.
One of the most interesting findings was the discovery of a strong GWAS signal that was mapped on a 1 Mb region of chromosome 3, which was associated with tolerance to leaf score injury under salt stress.The peak of this signal corresponded to Vigun03g144700.1, which encodes for a potassium channel.This potassium channel has been described to be activated upon salt stress in cowpea in order to enhance the transport of K + under salt stress in cowpea [28].Previous investigations have shown that salt-tolerant cowpea had a higher K + /Na + ratio in leaves [16,29].Therefore, the GWAS approach we used in this study has successfully targeted a gene that is involved in salt tolerance in cowpea.In addition, our previous research revealed a K + channel protein being involved in salt tolerance in a MAGIC population, which is in agreement with the GWAS result in this study.Potassium channel proteins have also been well-described for enhancing tolerance to salinity in other species.K + channel-related genes were shown to be upregulated under salt stress in tomato and soybean [30].
Genes encoding for NAD-dependent dehydratase have also been found in the vicinity of the significant SNPs associated with salt tolerance.These genes have been demonstrated to regulate stress in rice [31].A gene encoding for auxin efflux carrier was also found within the GWAS peaks.Auxin efflux proteins were reported to have a significant role in assisting Arabidopsis thaliana with regulating salt stress [32].The auxin efflux carriers regulate the variation in auxin flow during salt stress and are also involved in regulating meristem size for plants under salt stress.Results also indicated the involvement of a chlorophyllase gene in salt tolerance.However, there is no report yet highlighting the role of chlorophyllase in salt tolerance.We would suggest that chlorophyllase is a salt-susceptible gene since it is involved in chlorophyll degradation [33].Genes involved in vacuolar iron transporters were also identified.Our previous investigation on salt tolerance identified these genes in a MAGIC population.These transporters are also involved in salt tolerance [34].In soybean, the Na + /H + antiporter gene, GmCHX1, has been well described in conferring salt tolerance [35].A simple BLAST search showed that an orthologue of this gene can be found on chromosome 7 of cowpea.However, no strong GWAS peak was found on this chromosome.We could assume that this gene might be associated with a rare allele so that our GWAS approach failed to identify it.
A large number of molecular markers that are associated with cowpea salt tolerance have been identified in this study.A SNP validation is required prior to using these markers in a breeding program for Marker-Assisted Selection (MAS).The results from this investigation also contributed to a better understanding of the genetics of salt tolerance mechanism in cowpea.The candidate genes that were relevant to salt tolerance will be validated in further studies.Conducting the salt tolerance under greenhouse conditions could be a limitation of this study.However, to date, greenhouse phenotyping remains the most affordable and accurate way to evaluate salt tolerance since a lot of uncontrolled factors can occur during field screening.Therefore, repeating the experiments under field conditions could be a major challenge.

Plant Materials and Phenotyping
A total of 331 cowpea genotypes were evaluated for salt tolerance.This association panel consisted of breeding lines from the University of Arkansas, Fayetteville, and the University of California, Riverside [36], and Plant Introductions (PIs) from the U.S. Department of Agriculture (USDA) Germplasm Resources Information Network (GRIN) cowpea accessions and provided by the USDA Plant Genetic Resources Conservation Unit at Griffin, GA, USA.
Phenotypic data were collected in our previous project [25].Salt tolerance evaluation was conducted under greenhouse conditions at the University of Arkansas, Fayetteville, with average day/light temperatures of 26 • C/21 • C and an average daylight length of 14 hours.Salt tolerance was performed using a previously described methodology [37].Briefly, a total of four cowpea seedlings were established in pots filled with 100 g Sunshine Natural & Organic (Agawam, MA, USA).The experiments were conducted using a randomized complete block design (RCBD) with 2 replications within each block.A total of 4 blocks were used.Each pot corresponded to one replication.For each genotype, one pot was subjected to salt treatment, whereas another one was irrigated with deionized water and used as a control.A total of 12 pots were established on a rectangular plastic tray to facilitate irrigation.Salt treatment (NaCl) began at the first trifoliate leaf stage (V1 stage) [38] and was conducted by supplying a solution of 200 mM NaCl to each rectangular plastic tray [37].Two-thirds of pot's height was fully soaked with either deionized water or salt solution during irrigation [37].A salt-tolerant genotype ('09-529') and a salt-susceptible genotype (PI255774) [37] were used as checks.Data measurements were previously described, and the phenotypic data used for this study are from our previous report [25].From our previous phenotyping efforts, we found leaf injury score ranged between 1.4 to 9, leaf SPAD chlorophyll under salt stress varied from 6.4 to 39.9, and relative tolerance index for chlorophyll content ranged between 16.7 to 121.0.These phenotypic data indicate that a greater variability in salt tolerance traits exist among the cowpea accessions used in this study.4.2.Genotyping 4.2.1.DNA Extraction, Library Preparation, and Whole-Genome Resequencing Young cowpea leaves were harvested from one plant, and all seeds used during the experiments were derived from that one plant.Genomic DNA was extracted from freezedried young cowpea leaves using the CTAB (hexadecyltrimethyl ammonium bromide) protocol [39].Leaves were ground using a Mixer Mill MM 400 ® (Haan, Germany).DNA buffer was added to each sample.The mixture DNA buffer sample was centrifuged at 13,000 rpm for 10 min.Proteins were denatured by adding a solution of 1 mL of chloroform-isoamyl alcohol (24:1) to each sample.The addition of 1 mL of isopropanol helped DNA precipitate.In order to optimize DNA precipitation, samples were stored at −20 • C overnight.DNA pellets were washed using 70% and 90% ethanol.Washed DNA pellets were air dried.RNA was removed by adding 3 µl of RNAse to each sample.DNA was kept in a solution of 200 µL of 0.1× TE.DNA was quantity using a NanoDrop 200c spectrophotometer (Thermo SCIENTIFIC, Wilmington, DE, USA) and quality-checked on a 1%-agarose gel with ethidium bromide stain.

SNP Calling, Mapping, and Filtering
Reads were aligned to the cowpea reference genome [41] using SOAPaligner/soap2 (http://soap.genomics.org.cn/(accessed on 17 January 2020)).SNP calling was conducted using SOAPsnp v 1.05 [42].Accessions with more than 20% missing data were removed.Triallelic SNPs and those with more than 20% missing data were also removed.SNPs with a heterozygosity greater than 20% were removed as well.The minor allele frequency (MAF) threshold was 5%.GWAS was conducted using filtered SNPs.

Genome-Wide Association Study (GWAS)
GWAS was performed using Bayesian Information and Linkage Disequilibrium Iteratively Nested Keyway (BLINK) model [43].BLINK has been demonstrated to be statistically more powerful than the previously developed models [44].SNP was significant when above the FDR-adjusted threshold and computed in R (p < 3 × 10 −8 ).BLINK model was built upon the Fixed and Random Model Circulating Probability Unification (FarmCPU) model.In FarmCPU, markers were assumed to be evenly distributed across the genome, which was not necessarily true.BLINK used the LD information to relax this assumption.In addition, FarmCPU could be computationally intensive due to the random model part of its algorithm.The random model was replaced by a fixed model in BLINK.The two fixed effect models in BLINK are described below.

FEM (1)
FEM (2) with y i being the vector phenotype, M i1 , M i2 b 2 , . .., M ik the genotypes of k pseudo QTNs that were initially empty and with effects b 1 , b 2 , . .., b k , respectively, M ij being the j th genetic marker of the i th sample, and e i being the residual having a distribution with mean zero and a variance σ 2 e .Overlapping SNP markers between different traits were visualized using a Venn diagram and designed using the online software program accessible at http://jvenn.toulouse.inra.fr/app/example.html(accessed on 10 March 2020).

Candidate Gene Search and Synteny Analysis
By taking into account the number of SNPs involved in this study, the genome size of cowpea, and the average length of a gene within the cowpea genome, we investigated the annotated genes within 10 kb genomic region flanking an SNP using Phytozome v.13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 2 March 2020)).We considered annotated genes that were involved in plant physiology and/or tolerance to abiotic stress.Functional annotations of each annotated gene were obtained using Phytozome v. 13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 2 March 2020)).Coding sequences of the annotated genes relevant to plant physiology and/or tolerance to abiotic stress were extracted.The extracted sequences were used as query to perform BLASTx (https://blast.ncbi.nlm.nih.gov/Blast.cgi(accessed on 2 March 2020)) in order to obtain the amino acid sequence.Protein homolog search in other legumes such as soybean, common bean, and Medicago truncatula Gaertn was performed using the amino acid sequence.Only hits with similarity greater than 90% were considered.The tertiary structure of the amino acid sequence was predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 2 March 2020)).

Conclusions
In this study, strong GWAS peaks associated with leaf SPAD chlorophyll under salt stress, relative tolerance index for chlorophyll, and tolerance to leaf injury score under salt stress were identified.A total of 65, 60, and 1667 significant SNPs were found to be associated with leaf SPAD chlorophyll under salt stress, relative tolerance index for chlorophyll, and tolerance to leaf injury score under salt stress, respectively.Leaf SPAD chlorophyll under salt stress was characterized by a strong candidate locus by an 890 kb region of chromosome 2. Two candidate loci were found to be associated with the relative tolerance index for chlorophyll and mapped on chromosomes 1 and 2. A strong candidate locus defined by a 1-Mb region of chromosome 3 was associated with tolerance to leaf injury score under salt stress in cowpea.The results from this study could be used in cowpea breeding through Marker-Assisted Selection (MAS).To the best of our knowledge, this is the first report on cowpea GWAS using whole genome resequencing data.

Figure 1 .
Figure 1.Manhattan plot for leaf SPAD chlorophyll under salt stress in cowpea.The solid black and grey dots represent the SNPs.The x-axis is the chromosome number, and the y-axis is the LOD or −log10 of the p-value.The horizontal red and blue bars are two different LOD thresholds.Below the Manhattan plot are gene IDs from phytozome v.13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020)) corresponding to the significant locus on chromosome 1.Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast.cgi(accessed on 3 March 2020)).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 3 March 2020)) and presented on the left-hand side in the above figure.

Figure 1 .
Figure 1.Manhattan plot for leaf SPAD chlorophyll under salt stress in cowpea.The solid black and grey dots represent the SNPs.The x-axis is the chromosome number, and the y-axis is the LOD or −log10 of the p-value.The horizontal red and blue bars are two different LOD thresholds.Below the Manhattan plot are gene IDs from phytozome v.13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020)) corresponding to the significant locus on chromosome 1.Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast.cgi(accessed on 3 March 2020)).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 3 March 2020)) and presented on the left-hand side in the above figure.

Figure 2 .Figure 2 .
Figure 2. Manhattan plot for leaf SPAD chlorophyll under salt stress in cowpea.The solid black and grey dots represent the SNPs.The x-axis is the chromosome number, and the y-axis is the LOD or −log10 of the p-value.The horizontal red and blue bars are two different LOD thresholds.Below the Manhattan plot are gene IDs from phytozome v.13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020)) corresponding to the significant locus on chromosome 2. Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast.cgi(accessed on 3 March 2020)).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 3 March 2020)) and presented on the left-hand side in the above figure.Figure 2. Manhattan plot for leaf SPAD chlorophyll under salt stress in cowpea.The solid black and grey dots represent the SNPs.The x-axis is the chromosome number, and the y-axis is the LOD or −log10 of the p-value.The horizontal red and blue bars are two different LOD thresholds.Below the Manhattan plot are gene IDs from phytozome v.13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020)) corresponding to the significant locus on chromosome 2. Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast.cgi(accessed on 3 March 2020)).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 3 March 2020)) and presented on the left-hand side in the above figure.

Figure 3 .
Figure 3. Manhattan plot relative tolerance index for leaf SPAD chlorophyll for salt stress in cowpea.The solid black and grey dots represent the SNPs.The x-axis is the chromosome number, and the yaxis is the LOD or −log10 of the p-value.The horizontal red and blue bars are two different LOD thresholds.Below the Manhattan plot are gene IDs from phytozome v.13

Figure 3 .
Figure 3. Manhattan plot relative tolerance index for leaf SPAD chlorophyll for salt stress in cowpea.The solid black and grey dots represent the SNPs.The x-axis is the chromosome number, and the y-axis is the LOD or −log10 of the p-value.The horizontal red and blue bars are two different LOD thresholds.Below the Manhattan plot are gene IDs from phytozome v.13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020)) corresponding to the significant locus on chromosome 1.Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast (accessed on 3 March 2020)).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 3 March 2020)) and presented on the left-hand side in the above figure.
(https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020)) corresponding to the significant locus on chromosome 1.Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast(accessed on 3 March 2020)).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 3 March 2020)) and presented on the left-hand side in the above figure.

Figure 4 .
Figure 4. Manhattan plot relative tolerance index for leaf SPAD chlorophyll for salt stress in cowpea.The solid black and grey dots represent the SNPs.The x-axis is the chromosome number, and the yaxis is the LOD or −log10 of the p-value.The horizontal red and blue bars are two different LOD thresholds.Below the Manhattan plot are gene IDs from phytozome v.13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020)) corresponding to the significant locus on chromosome 2. Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast.cgi(accessed on 3 March 2020)).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 3 March 2020)) and presented on the left-hand side in the above figure.

Figure 4 .
Figure 4. Manhattan plot relative tolerance index for leaf SPAD chlorophyll for salt stress in cowpea.The solid black and grey dots represent the SNPs.The x-axis is the chromosome number, and the y-axis is the LOD or −log10 of the p-value.The horizontal red and blue bars are two different LOD thresholds.Below the Manhattan plot are gene IDs from phytozome v.13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020)) corresponding to the significant locus on chromosome 2. Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast.cgi(accessed on 3 March 2020)).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 3 March 2020)) and presented on the left-hand side in the above figure.

Figure 5 .
Figure 5. Manhattan plot relative tolerance index for leaf SPAD chlorophyll for salt stress in cowpea.The solid black and grey dots represent the SNPs.The x-axis is the chromosome number, and the yaxis is the LOD or −log10 of the p-value.The horizontal red and blue bars are two different LOD thresholds.Below the Manhattan plot are gene IDs from phytozome v.13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020)) corresponding to the significant locus on chromosomes 3 and 4. Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast.cgi(accessed on 3 March 2020)).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 3 March 2020)) and presented on the left-hand side in the above figure.

Figure 5 .
Figure 5. Manhattan plot relative tolerance index for leaf SPAD chlorophyll for salt stress in cowpea.The solid black and grey dots represent the SNPs.The x-axis is the chromosome number, and the y-axis is the LOD or −log10 of the p-value.The horizontal red and blue bars are two different LOD thresholds.Below the Manhattan plot are gene IDs from phytozome v.13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020)) corresponding to the significant locus on chromosomes 3 and 4. Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast.cgi(accessed on 3 March 2020)).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 3 March 2020)) and presented on the left-hand side in the above figure.
zome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020)) corresponding to the significant locus on chromosomes 3 and 4. Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast.cgi(accessed on 3 March 2020)).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 3 March 2020)) and presented on the left-hand side in the above figure.

Figure 6 .
Figure 6.Manhattan plot relative tolerance index for leaf SPAD chlorophyll for salt stress in cowpea.The solid black and grey dots represent the SNPs.The x-axis is the chromosome number, and the y-axis is the LOD or −log10 of the p-value.The horizontal red and blue bars are two different LOD thresholds.Below the Manhattan plot are gene IDs from phytozome v.13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020)) corresponding to the significant locus on chromosome 8.Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast.cgi(accessed on 3 March 2020)).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 3 March 2020)) and presented on the left-hand side in the above figure.

Figure 7 .
Figure 7. Manhattan plot for tolerance to leaf injury score under salt stress in cowpea.The solid black and grey dots represent the SNPs.The x-axis is the chromosome number, and the y-axis is the LOD or −log10 of the p-value.The horizontal red and blue bars are two different LOD thresholds.Below the Manhattan plot are gene IDs from phytozome v.13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020)) corresponding to the significant locus on chromosome 1.Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast.cgi(accessed on 3 March 2020)).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 3 March 2020)) and presented on the left-hand side in the above figure.

Figure 7 .
Figure 7. Manhattan plot for tolerance to leaf injury score under salt stress in cowpea.The solid black and grey dots represent the SNPs.The x-axis is the chromosome number, and the y-axis is the LOD or −log10 of the p-value.The horizontal red and blue bars are two different LOD thresholds.Below the Manhattan plot are gene IDs from phytozome v.13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020)) corresponding to the significant locus on chromosome 1.Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast.cgi(accessed on 3 March 2020)).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 3 March 2020)) and presented on the left-hand side in the above figure.

Figure 7 .
Figure 7. Manhattan plot for tolerance to leaf injury score under salt stress in cowpea.The solid black and grey dots represent the SNPs.The x-axis is the chromosome number, and the y-axis is the LOD or −log10 of the p-value.The horizontal red and blue bars are two different LOD thresholds.Below the Manhattan plot are gene IDs from phytozome v.13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020) corresponding to the significant locus on chromosome 1.Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast.cgiaccessed on 3 March 2020).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/accessed on 3 March 2020) and presented on the left-hand side in the above figure.

Figure 8 .
Figure 8. Manhattan plot for tolerance to leaf injury score under salt stress in cowpea.The solid black and grey dots represent the SNPs.The x-axis is the chromosome number, and the y-axis is the LOD or −log10 of the p-value.The horizontal red and blue bars are two different LOD thresholds.Below the Manhattan plot are gene IDs from phytozome v.13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020)) corresponding to the significant locus on chromosome 2. Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast.cgi(accessed on 3 March 2020)).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 3 March 2020)) and presented on the left-hand side in the above figure.

Figure 9 .
Figure 9. Manhattan plot for tolerance to leaf injury score under salt stress in cowpea.The solid black and grey dots represent the SNPs.The x-axis is the chromosome number, and the y-axis is the LOD or −log10 of the p-value.The horizontal red and blue bars are two different LOD thresholds.Below the Manhattan plot are gene IDs from phytozome v.13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020)) corresponding to the significant locus on chromosome 3. Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast.cgi(accessed on 3 March 2020)).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 3 March 2020)) and presented on the left-hand side in the above figure.

Figure 9 .
Figure 9. Manhattan plot for tolerance to leaf injury score under salt stress in cowpea.The solid black and grey dots represent the SNPs.The x-axis is the chromosome number, and the y-axis is the LOD or −log10 of the p-value.The horizontal red and blue bars are two different LOD thresholds.Below the Manhattan plot are gene IDs from phytozome v.13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020)) corresponding to the significant locus on chromosome 3. Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast.cgi(accessed on 3 March 2020)).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 3 March 2020)) and presented on the left-hand side in the above figure.

Figure 10 .
Figure 10.Manhattan plot for tolerance to leaf injury score under salt stress in cowpea.The solid black and grey dots represent the SNPs.The x-axis is the chromosome number, and the y-axis is the LOD or −log10 of the p-value.The horizontal red and blue bars are two different LOD thresholds.Below the Manhattan plot are gene IDs from phytozome v.13 (https://phytozome.jgi.doe.gov/pz/portal.html(accessed on 3 March 2020)) corresponding to the significant locus on chromosome 4. Codifying sequences of the gene IDs whose functions were related to drought stress were extracted and converted to amino acid sequences using BLASTX (https://blast.ncbi.nlm.nih.gov/Blast.cgi(accessed on 3 March 2020)).Tertiary structures of the proteins/polypeptides derived from BLASTX were predicted using SWISS-MODEL (https://swissmodel.expasy.org/(accessed on 3 March 2020)) and presented on the left-hand side in the above figure.

21 Figure 11 .
Figure 11.(A) Venn diagram showing the number of overlapping significant SNPs associated with leaf SPAD chlorophyll under salt stress (S_Chloro), relative tolerance index for leaf SPAD chlorophyll (S_RTI), and leaf injury score under salt stress (S_Score) in cowpea.(B) Venn diagram showing the number of unique functional annotations for candidate genes associated with leaf SPAD chlorophyll under salt stress (S_Chloro), relative tolerance index for leaf SPAD chlorophyll (S_RTI), and leaf injury score under salt stress (S_Score) in cowpea.Venn diagrams were established using http://jvenn.toulouse.inra.fr/app/example.html(accessed on 10 March 2020).

Figure 11 .
Figure 11.(A) Venn diagram showing the number of overlapping significant SNPs associated with leaf SPAD chlorophyll under salt stress (S_Chloro), relative tolerance index for leaf SPAD chlorophyll (S_RTI), and leaf injury score under salt stress (S_Score) in cowpea.(B) Venn diagram showing the number of unique functional annotations for candidate genes associated with leaf SPAD chlorophyll under salt stress (S_Chloro), relative tolerance index for leaf SPAD chlorophyll (S_RTI), and leaf injury score under salt stress (S_Score) in cowpea.Venn diagrams were established using http://jvenn.toulouse.inra.fr/app/example.html(accessed on 10 March 2020).

Table 1 .
List of significant SNPs close to candidate genes and associated with leaf SPAD chlorophyll under salt stress, relative tolerance index for chlorophyll, and leaf injury score under salt stress in cowpea.SNP, CHR, BP, Pval, and LOD refers to SNP_ID, chromosome number, physical location (in bp), p-value, and -log10 of p-value (LOD), respectively.Gene_ID and functional annotations were obtained from Pythozome v.13.