Genetic architecture of salt tolerance in a Multi-Parent Advanced Generation Inter-Cross (MAGIC) cowpea population

Background Previous reports have shown that soil salinity is a growing threat to cowpea production, and thus the need for breeding salt-tolerant cowpea cultivars. A total of 234 Multi-Parent Advanced Generation Inter-Cross (MAGIC) lines along with their 8 founders were evaluated for salt tolerance under greenhouse conditions. The objectives of this study were to evaluate salt tolerance in a multi-parent advanced generation inter-cross (MAGIC) cowpea population, to identify single nucleotide polymorphism (SNP) markers associated with salt tolerance, and to assess the accuracy of genomic selection (GS) in predicting salt tolerance, and to explore possible epistatic interactions affecting salt tolerance in cowpea. Phenotyping was validated through the use of salt-tolerant and salt-susceptible controls that were previously reported. Genome-wide association study (GWAS) was conducted using a total of 32,047 filtered SNPs. The epistatic interaction analysis was conducted using the PLINK platform. Results Results indicated that: (1) large variation in traits evaluated for salt tolerance was identified among the MAGIC lines, (2) a total of 7, 2, 18, 18, 3, 2, 5, 1, and 23 were associated with number of dead plants, salt injury score, leaf SPAD chlorophyll under salt treatment, relative tolerance index for leaf SPAD chlorophyll, fresh leaf biomass under salt treatment, relative tolerance index for fresh leaf biomass, relative tolerance index for fresh stem biomass, relative tolerance index for the total above-ground fresh biomass, and relative tolerance index for plant height, respectively, with overlapping SNP markers between traits, (3) candidate genes encoding for proteins involved in ion transport such as Na+/Ca2+ K+ independent exchanger and H+/oligopeptide symporter were identified, and (4) epistatic interactions were identified. Conclusions These results will have direct applications in breeding programs aiming at improving salt tolerance in cowpea through marker-assisted selection. To the best of our knowledge, this study was one of the earliest reports using a MAGIC population to investigate the genetic architecture of salt tolerance in cowpea. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-022-08332-y.

million hectares of croplands [3]. Cowpea is a legume that has a multipurpose use. It provides an excellent and affordable source of protein to human [4]. Cowpea seeds contain nutrients that are necessary to human's heath. One hundred g of cowpea seed contains 6.8 mg of iron, 4.1 mg of zinc, 1.5 mg of manganese, 510 mg of phosphorus, and 1430 mg of potassium [5]. The significant amount of antioxidant compounds within cowpea seeds provides additional nutritional value that would be of interest when incorporated into the diet [6,7]. In addition to significantly contributing to enhancing the human's diet, cowpea hay could be used to supplement low quality feed for livestock, a prevalent practice in sub-Sahara Africa [2]. Cowpea also contributes to effective ecosystem management by limiting soil erosion. In fact, with its excellent root architecture plus nitrogen fixation, cowpea can be used as a cover crop that has attracted considerable attention in recent years [8].
Despite of the aforementioned benefits, cowpea production can be substantially limited by abiotic stresses also including soil salinity. Salinity has been reported to increasingly affecting agricultural production worldwide and contributing to an annual loss of 12 billion US dollars [9,10]. Soil salinity has resulted from the accumulation of cations consisting of K + , Mg 2+ , Ca 2+ , and Na + and anions such as NO 3 − , HCO 3 − , SO 4 2− , and Cl − within the soil profile [11]. Soil salinity affects more than 19.6 million ha of croplands in the U.S. and areas facing salinity-related issues have increased [12]. Cowpea cultivation is common in semi-arid areas since cowpea has a better capability to withstand a limited water condition [13]. However, earlier reports suggested that the limited rainfall occurring in semi-arid areas significantly contributed to the salt-related compounds not being effectively leached out from the soil profile, which can exacerbate the effects of salinity on cowpea grown in semi-arid regions [14]. Salinity is also increased by the use of poor quality irrigation water. In the U.S., cowpea cultivation is prevalent in the southern regions [15]. However, irrigation from groundwater in the southern U.S. accounts for more than 66% of the water source used for agricultural activities and can contain up to 1639 mg of Cl − per L of water [16,17]. A sodium chloride (NaCl) concentration greater than 90 mM, releasing around 526 mg/L of Cl − , could significantly reduce cowpea yield [18]. Therefore, salinity could limit cowpea production in relevant southern U.S. areas. Significant cowpea production can also be found in western U.S. in addition to the increasing interest in the use of cowpea as cover crop in this part of the country [8]. The Coachella Valley of California has also been increasingly impacted by salinity, limiting cowpea cultivation expansion in western U.S. [8,19]. Salinity can also be increased by the overuse of fertilizers or natural factors such as rock weathering [20].
Salinity affects most of development and growth stages of cowpea with germination and seedling stages being the most sensitive stages [21,22]. Salinity can completely suppress cowpea germination and lead to plant death in cowpea seedlings [22]. In addition, high salt ion concentrations will result in significant reductions in height, biomass, and chlorophyll content, leading to serious physiological impairment within cowpea plants [21]. Breeding for salt-tolerant cowpea cultivars would be one of the most affordable ways to limit the negative effects of salinity on cowpea cultivation. Significant efforts towards investigating salt tolerance in cowpea have been conducted in relatively recent years. Salt tolerance at the germination stage among 151 diverse cowpea accessions have been reported [22]. Another study have identified promising cowpea lines that better withstand salt stress at seedling stage [21]. Molecular markers have substantially assisted plant breeders with rapidly developing cultivars. Our previous article reported the first molecular markers associated with salt tolerance in cowpea [23]. Three SNPs (Scaffold87490_622, Scaffold87490_630, and C35017374_128) were found to be associated with salt tolerance at both seedling and germination stages, and other seven SNPs (Scaffold93827_270, Scaf-fold68489_600, Scaffold87490_633, Scaffold87490_640, Scaffold82042_3387, C35069468_1916, and Scaf-fold93942_1089) were reported to be associated with seedling stage-specific salt tolerance in cowpea [23]. The aforementioned research was carried out on an association panel consisting of a diverse cowpea germplasm but having a limited population size, which would reduce the likelihood of finding rare alleles that potentially affect salt tolerance. This can be addressed by conducting a genome-wide association analysis (GWAS) on a multi-parent advanced generation inter-cross (MAGIC) population, as it can increase the frequency of rare alleles while providing significant recombination between the chromosomal Sects. [24][25][26][27]. In this study, we used the Bayesian Information and Linkage Disequilibrium Iteratively Nested Keyway (BLINK) model to conduct GWAS. BLINK was built upon the Fixed and Random Model Circulating Probability Unification (FarmCPU) model. In FarmCPU, markers are assumed to be evenly distributed across the genome. However, such an assumption could be easily violated. BLINK relaxed this assumption by incorporating the LD information. The random effect model (REM) part in FarmCPU, which was computationally heavy, was replaced by a second fixed effect model (FEM) in BLINK.
The first MAGIC cowpea population [28] has been developed by using eight founder parents having desirable agronomic traits such as high yield, drought tolerance, and resistance to diseases and insects [28]. However, salt tolerance was not investigated for this MAGIC population while salinity has become an increasing threat to cowpea production worldwide. In addition, genomic selection for salt tolerance remains very limited in cowpea improvement. Investigating the epistatic interactions between markers will further advance our understanding of the genetic mechanism affecting salt tolerance in cowpea. To date, no studies have explored the possibility of epistatic interactions governing salt tolerance in cowpea. The recent advances made in quantitative genetics and modelling permit the identification of epistatic interactions across the plant genome. The study of epistatic interaction between markers can also provide insights to possible gene pathways affecting salt tolerance. These pathways have been demonstrated to affect salt tolerance. For example, in rice, genes involved in membrane sensor proteins, signaling pathways, and electron carriers were found to confer salt tolerance [29]. In cotton, genes involved in the ABA signaling pathway have been reported to enhance salt tolerance [30]. Other reports also highlighted that the salt overly sensitive (SOS) is a key pathway in regulating salt stress in plants [31]. Genes involved in the jasmonate pathway have also been shown to contribute to salt tolerance [32]. The objectives of this study were to evaluate salt tolerance in the cowpea MAGIC population, to identify SNP markers associated with this trait, to assess the accuracy of GS for salt tolerance, and to explore the potential epistatic interactions affecting salt tolerance between SNP markers.

Phenotypic data
The average number of dead plants per pot varied from 0.0 to 3.0, with an average of 1.0 and a standard deviation of 0.7 (Table S1). The distribution of the average number of dead plants per pot was right-skewed ( Fig. 1 A). A significant difference in the average number of dead plants per pot was found among the lines (F-value=15.3, p-value<0.0001) ( Table 1), which was expected. Despite the significant line X block interaction effect, the main factor (line) was still analyzed since assessing salt tolerance between lines was the main purpose of the phenotypic evaluation in this study. Of the 242 genotypes evaluated for salt tolerance, 45 did not have any dead plants across 4 replications. In addition, variation in the average number of dead plants per pot was identified as shown in Fig. 1 A. Interestingly, none of the cowpea parents were among the worst 45 lines with plant death. The  The 8  founders were P1: CB27, P2: IT00K-1263, P3: IT82E-18, P4: IT84S-2049, P5: IT84S-2246, P6: IT89KD-288, P7: IT93K-503-1, and P8: Suvita 2  Table 1 ANOVA for the MAGIC population evaluated under salt tolerance. The evaluated traits were the average number of dead plants per pot (Dead_plants), leaf injury score (Salt_score), SPAD chlorophyll under salt treatment (S_Chloro), relative tolerance index for SPAD chlorophyll (RTI_C), fresh leaf biomass under salt treatment (S_Leaf ), relative tolerance index for fresh leaf biomass (RTI_FL), relative tolerance index for fresh stem biomass (RTI_FS), relative tolerance index for total above-ground fresh biomass (RTI_FB), and relative tolerance index for plant height (RTI_H)  Table 2). The broad sense heritability for the average number of dead plants per pot was 74.2%. Leaf injury score was approximately normally distributed (Fig. 1B) and ranged between 0.5 and 6.5, with an average of 3.6 and a standard deviation of 1.1 ( Table S1). The lines differed significantly in leaf injury (F-value=13.2, p-value<0.0001) ( Table 1). Lines with the lowest leaf injury score were MAGIC208 (0.  (Table 2),.Lines with the highest leaf injury score were MAGIC259 (6.0), MAGIC298 (6.0), MAGIC194 (6.0), MAGIC048 (6.3), MAGIC092 (6.5) ( Table 2), which were the most susceptible genotypes in terms of leaf injury score. None of the MAGIC parents were among the most tolerant and the most susceptible groups. The parent that was the most tolerant to salt stress in terms of leaf injury score were IT00K-1263 (3.3), whereas the one that was the most susceptible was Table 2 List of top 5 genotypes and 5 least performers for average number of dead plants per plot (DeadPlants), the leaf injury score under salt treatment (Score), SPAD chlorophyll under salt treatment (StressSPADChloro), relative tolerance index for SPAD chlorophyll (RTI_C), fresh leaf biomass under salt treatment (StressLeaf ), relative tolerance index for fresh leaf biomass (RTI_FL), relative tolerance index for fresh stem biomass (RTI_FS), relative tolerance index for total above-ground fresh biomass (RTI_FB), and relative tolerance index for plant height (RTI_H). Sd represents the standard deviation across 4 replications IT89KD-288 (5.8) (Fig. 1B). The broad sense heritability (H) for leaf injury score under salt stress was 72.6%. The distribution of leaf SPAD chlorophyll under salt treatment was approximately normally distributed as shown in Fig. 1 C. Leaf SPAD chlorophyll was significantly different among the MAGIC lines (F-value=45.2, p-value<0.0001) ( Table 1). The lines with the highest leaf SPAD chlorophyll under salt stress were MAGIC208 (51.  Table 2). The MAGIC parent with the highest leaf SPAD chlorophyll under salt treatment was IT84S-2246 (13.6), which was the most susceptible parent in terms of leaf SPAD chlorophyll under salt treatment. The MAGIC parent with the highest leaf SPAD chlorophyll under salt stress was IT93K-503-1 (21.9) ( Table 2). The broad sense heritability (H) for leaf SPAD chlorophyll under salt treatment was 78.9%.
Fresh leaf biomass under salt stress is also a good character for assessing salt tolerance in cowpea at the seedling stage. In this study, fresh leaf biomass of cowpea plants under salt treatment was approximately normally distributed (Fig. 1E). Fresh leaf biomass ranged between 0.5 and 4.2 g, with an average of 2.1 g and a standard deviation of 0.7 g (Table S2). Under salt stress, a significant difference in fresh leaf biomass was identified among the lines (F-value=11.9, p-value<0.0001) ( Table 1). Lines with the highest fresh leaf biomass under salt stress were MAGIC208 (4.2 g), MAGIC336 (3.8 g), MAGIC271 (3.8 g), MAGIC187 (3.8 g), and MAGIC027 (3.8 g), whereas those with the lowest fresh leaf biomass under salt stress were Suvita 2 (0.7 g), MAGIC073 (0.6 g), MAGIC048 (0.6 g), MAGIC092 (0.5 g), and IT84S_2246 (0.5 g) ( Table 2). The broad sense (H) heritability for fresh leaf biomass of cowpea plants grown under salt treatment was 61.3%.
The distribution of relative tolerance index for plant height (RTI_H) was approximately normal (Fig. 1I). RTI_H varied from 54.6 to 89.5%, with an average of 73.1% and a standard deviation of 6.0% (Table S3). A significant difference in RTI_H was identified among the lines (F-value=6.9, p-value<0.0001) ( Table 1

Correlation analysis
The average number of dead plants per pot was strongly correlated with the salt injury score (r=0.9, P < 0.001). In addition, a strong and negative correlation was found between the average number of dead plants per pot and the leaf SPAD chlorophyll under salt stress (r=-0.8, P < 0.001), and between the average number of dead plants per pot and the relative tolerance index for leaf SPAD chlorophyll (r=-0.8, P < 0.001) ( Table 3). A high and negative correlation was also identified between the average number of dead plants per pot and fresh leaf biomass under salt stress (r=-0.6, P < 0.001). The average number of dead plants per pot was significantly correlated with the relative tolerance index for fresh stem biomass (r=-0.20, P < 0.001), relative tolerance index for total above-ground fresh biomass (r=-0.30, 0.001), and relative tolerance index for plant height (r=-0.10, nonsignificant) ( Table 3). The relative tolerance index for leaf SPAD chlorophyll was moderately correlated with fresh leaf biomass under salt stress (r=0.50, P < 0.001), but the relative tolerance index for leaf SPAD chlorophyll was not correlated with the relative tolerance index for plant height (r=-0.10, non-significant), indicating that the mechanism for tolerance to plant height reduction and leaf SPAD chlorophyll reduction under salt stress could be different. The trait having the highest correlation with relative tolerance index for plant height was fresh stem biomass (r=0.40, P < 0.001) ( Table 3).
The pairwise relationship that was based on the Pearson's correlation coefficient for the traits evaluated under salt stress was visualized used a chord diagram (Fig. 2). The thicker the link between traits the higher the Person's correlation coefficient. Traits with the thickest link end were the average number of dead plants per pot, leaf injury score, leaf SPAD chlorophyll under salt stress, and relative tolerance index for leaf SPAD chlorophyll (Fig. 2), suggesting the possibility of common pathway(s) for salt tolerance mechanism for these traits. Traits with the thinnest link end were relative tolerance index for fresh stem and plant height, indicating that the mechanism for salt tolerance could be independent from the other traits evaluated in this study.

GWAS and candidate gene identification
A total of seven significant SNPs were identified to be associated with the average number of dead plants per pot (Table 4), among which three SNPs were located on chromosome 3 and four SNPs on chromosome 7 (Fig. 3). The Table 3 Persons' correlation coefficients for the traits evaluated for salt tolerance in a MAGIC population. Traits consisted of average number of dead plants per plot (DeadPlants), the leaf injury score under salt treatment (Score), SPAD chlorophyll under salt treatment (StressSPADChloro), relative tolerance index for SPAD chlorophyll (RTI_C), fresh leaf biomass under salt treatment (StressLeaf ), relative tolerance index for fresh leaf biomass (RTI_FL), relative tolerance index for fresh stem biomass (RTI_FS), relative tolerance index for total above-ground fresh biomass (RTI_FB), and relative tolerance index (RTI) for plant height (RTI_Height) *, **, *** indicates a significant correlation at P < 0.05, 0.01, and 0.001, respectively

RTI_FL
GWAS for RTI_C identified a total of 17 significant SNPs (Table 4). These SNPs were the ones that were associated with leaf SPAD chlorophyll under salt stress and were located within the 4.2-Mb genomic region of chromosome 3 (Fig. 3), suggesting a high likelihood Fig. 3 Manhattan plots for genome-wide association study (GWAS) corresponding to the average number of dead plants per pot (Dead), leaf injury score (Score), leaf SPAD chlorophyll under salt stress (S_Chloro), relative tolerance index for leaf SPAD chlorophyll (RTI_C), fresh leaf biomass under salt stress (S_Leaf ), relative tolerance index for fresh leaf biomass (RTI_FL), relative tolerance for fresh stem biomass (RTI_FS), relative tolerance index for total above-ground fresh biomass (RTI_FB), and relative tolerance index for plant height (RTI_H). For each Manhattan plot, the x-axis represents the chromosome number and the y-axis indicates the -log 10 (p) where is the p-value corresponding to each SNP after running BLINK. The horizontal red line indicates the LOD threshold Table 4 List of SNPs significantly associated with the traits evaluated under salt tolerance in a MAGIC cowpea population, chromosome and physical position (bp) of each SNP, LOD (-log10(p-value)), minor allele frequency MAF (%), annotated genes found within the 20-kb region flanking each significant SNP, and functional annotations for each gene ID. LOD threshold was greater or equal to 3.5. If no SNPs were above the threshold, the top 3 SNPs with the highest LOD were listed in below  Two SNPs were found to be significantly associated with RTI_FL (Fig. 3). These SNPs were also identified to be associated with fresh leaf biomass under salt stress, RTI_C, and leaf SPAD chlorophyll under salt stress (Table 4). While no SNPs were above the chosen LOD threshold  One SNP, 2_33574, was significantly associated with RTI_FB (Fig. 3). This SNP was located at 579,544 bp on chromosome 5 ( Table 4). The annotated genes found in the vicinity of this SNP were Vigun05g006800.1, Vigun05g006700.1, Vigun05g006600.1, and Vigun05g006500.1. No functional annotations were found for Vigun05g006600.1. Functional annotations for Vigun05g006800.1, Vigun05g006700.1, and Vigun05g006500.1 were Mannose-6-phosphate isomerase, alpha/beta hydrolase fold-containing protein, and neoxanthin biosynthesis, respectively. GWAS suggested a strong candidate locus associated with RTI_H) (Fig. 3). This genomic region harbored a total of 23 significant SNPs and were mapped on a 3.4-Mb region of chromosome 3 (  (Table 4). Functional annotations associated with the candidate genes were O-methyltransferase-related, protein transport protein Sect. 23, peptidyl-prolyl cis-trans isomerase, cystatin-C, phospholipases, dolichol-phosphate mannosyltransferase, IQ-domain 9 protein, mutt-nudix-related, magnesium chelatase subunit I, ionotropic glutamate receptor, apoptosis inhibitor 5, peroxidase 19, triacylglycerol degradation, cytochrome P450, microfibrilassociated protein, suberin monomers biosynthesis, homoserine dehydrogenase, and beta-galactosidase 9 ( Table 4).

Overlapping SNPs between traits
Common SNP markers were identified between the traits evaluated for salt tolerance in this MAGIC cowpea population. Out of the SNP markers associated with the average number of dead plants per pot, a total of 3 were found to be associated with both leaf SPAD chlorophyll under salt treatment (S_Chloro) and relative tolerance index for leaf SPAD chlorophyll (RTI_C) (Fig. 4), indicating that there could be a common pathway for salt tolerance based on these traits. A total of 14 significant SNPs were overlapping between S_Chloro and RTI_C (Fig. 4). Interestingly, none of the significant SNP markers associated with relative tolerance index for plant height (RTI_H) overlapped with any SNP markers associated with other traits (Fig. 4), suggesting that the mechanism for salt tolerance based on RTI_H could be independent. Using a Venn diagram with more than 5 sets would be difficult to visualize, so the Venn diagram (Fig. 4) did not include the data for leaf injury score (Score), relative tolerance index for fresh stem biomass (RTI_FS), relative tolerance index for the total above-ground fresh biomass (RTI_FB), fresh leaf biomass under salt stress (S_Leaf ), and relative tolerance index for fresh leaf biomass (RTI_ FL). None of the SNP makers associated with Score overlapped with any SNP makers associated with other traits (Table 4). Similar results were found for RTI_FS and RTI_FB. One SNP associated with S_Leaf, 2_28348, overlapped with one SNP associated with RTI_FL, S_Chloro, Dead, and RTI_C ( Table 4). The SNP 2_27478, associated with S_Leaf, was also associated with RTI_C, S_Chloro, and RTI_FL (Table 4). These results indicated that there could be a common pathway for salt tolerance between S_Leaf, RTI_FL, S_Chloro, Dead, and RTI_C.
All pairwise epistatic interactions found for the average number of dead plants per pot were between chromosomes with chromosomes 9 and 11 having the highest number of significant epistasis (Fig. 5 A). However, these epistasis-rich regions had SNPs with low LOD values. The genomic region of chromosome 7 that harbors some of the significant SNP markers associated with the average number of dead plants per pot was in epistasis with some SNPs found at the beginning of chromosome 8 (Fig. 5 A). No significant interaction was identified between SNPs located within the two candidate loci (one on chromosome 3 and the other on chromosome 7) associated with the average number of dead plants per pot. Similar results were found for leaf injury score where no epistatic interactions were identified between the significant SNPs associated this trait (Fig. 5B). The chromosomes with the highest number of epistatic interactions were chromosome 3 and chromosome 8 (Fig. 5B). Interestingly, most of significant epistatic interactions for leaf injury score appeared to be located towards both ends of the chromosome as shown in Fig. 5B. The epistasis analysis results for S_Chloro were particularly complex in which significant SNP markers associated with this trait, Fig. 4 Venn diagram showing the overlapping significant SNP markers between the average number of dead plants per pot (Dead), leaf SPAD chlorophyll under salt treatment (S_C), relative tolerance index for leaf SPAD chlorophyll (RTI_Chloro), and relative tolerance index for plant height (RTI_H). Venn diagrams were established using the online software program that is accessible at http:// jvenn. toulo use. inra. fr/ app/ examp le. html which were located on chromosome 3, were in epistatic interaction with SNPs located on chromosomes 2, 8, and 11 ( Fig. 5 C).
Results indicated a within-chromosome epistatic interaction (chromosome 4) for RTI_C (Fig. 5D). The pattern of epistasis for RTI_C was very similar to that of S_ Chloro (Fig. 5 C and D), which was expected since these traits were highly correlated ( Table 3). The chromosomes with the highest number of significant epistasis for S_ Leaf were 3 and 4 (Fig. 5E). None of the significant SNP markers associated with S_Leaf were in epistasis with any other SNPs. For RTI_FL, chromosomes 6 and 7 had the highest number of significant epistatic interactions and a within-chromosome epistasis was found on chromosome 7 (Fig. 6 A). The significant SNP markers associated with RTI_FL and found on chromosome 3 were not in epistatic interaction with any other SNPs (Fig. 6 A).
Epistatic interactions were also identified for RTI_FS. The chromosomes with the highest epistatic interaction were 1, 6, and 11 (Fig. 6B). The significant SNPs associated with RTI_FS and located on chromosome 4 were in epistatic interaction with some low-LOD SNPs on chromosome 7 (Fig. 6B). However, the other candidate locus containing significant SNPs and located on chromosome 7 was not in epistatic interaction with any genomic regions. RTI_FB had the highest number of significant epistatic interactions among all traits evaluated for salt tolerance in this study. The chromosomes with the highest number of significant epistatic interactions for RTI_FB were 3, 6, and 10 ( Fig. 6 C). The significant SNP markers associated with RTI_FB and mapped on chromosome 5 were in epistatic interaction with some low-LOD SNPs on chromosome 6. No within-chromosome epistatic interactions were identified for RTI_FB. The significant SNP markers associated with RTI_H and located on chromosome 3 were not in epistatic interaction with any other SNPs as shown in Fig. 6D. One within-chromosome epistatic interaction was found on chromosome 7.

Discussions
The large phenotypic variation and heritabilities observed in 234 MAGIC lines and 8 founder parents suggest that this population would be sufficient for salt tolerance genetic studies. To the best of our knowledge, this study Fig. 6 Circos plots showing the significant pairwise epistatic interactions between SNPs. On each circos plot, the outermost layer represents the 11 chromosomes of cowpea and the length of each segment is proportional to the length of each chromosome. The innermost layer displays the SNPs used for conducting GWAS and each black dot represents one SNP. The width of the innermost layer is proportional to the LOD values of each SNP. The further from the center the black dot is, the higher the LOD is. Links within each circos plot show a significant epistatic interaction between two SNPs. Since the resolution of the chromosomal length is in Mb (outermost layer), two closely located pairs of pairwise epistatic interactions can be cofounded in the above figure, so the number of links might not reflect the actual number of pairwise epistatic interactions. A Relative tolerance index for fresh leaf biomass, B relative tolerance index for fresh stem biomass, C relative tolerance index for total above-ground fresh biomass, and D relative tolerance index for plant height was one of the earliest reports investigating salt tolerance based on a MAGIC population in cowpea. In addition, the population size for this study was larger than those used in previous reports investigating cowpea salt tolerance [21,23].
GWAS identified significant SNP markers associated with various traits evaluated under salt stress in this MAGIC cowpea population, expanding the list of SNP markers associated with other important traits that were previously reported for cowpea [7,[33][34][35]. The earliest study for salt tolerance also identified several SNPs associated with salt tolerance in cowpea, including Scaffold87490_622, Scaf-fold87490_630, C35017374_128, Scaffold93827_270, Scaffold68489_600, Scaffold87490_633, Scaf-fold87490_640, Scaffold82042_3387, C35069468_1916, and Scaffold93942_1089 [23]. These SNPs were identified by conducting GWAS based on a diversity panel of 155 cowpea lines and 1049 SNPs that were postulated from genotyping-by-sequencing. In the present study, we performed GWAS based on a larger panel genotyped with a large number of SNPs. In addition, most of the SNP markers identified in this study were within or in the vicinity of annotated genes whose functional annotations involved salt tolerance mechanisms, which provides confidence to our findings.
Various candidate genes encoding for proteins with possible functions in salt tolerance mechanisms have been identified. The candidate gene (Vigun01g093100.1) revealed a relationship between Na + /Ca2 + K + independent exchanger and salt tolerance in cowpea. The involvement of Na + /Ca2 + K + independent exchanger in salt tolerance has been well described in other species such as tomato and soybean [36]. Therefore, the SNP markers (2_13484 and 2_13485) found in the vicinity of this gene could be reliably used for screening salt tolerance in cowpea. In addition, a H+/oligopeptide symporter has been shown to be associated with Cl-dynamic under salt stress in soybean [37]. These results suggested that a common salt tolerance mechanism pathway could exist between soybean and cowpea. Calcium-dependent protein kinases have also been identified to be associated with salt tolerance based on our data. Another study showed that calciumdependent protein kinases are important in regulating responses to salt stress in cotton [38]. These proteins play a role in stress signaling. Transcripts encoding for calcium-dependent protein kinases were also reported to be induced at early stage of salt stress in cotton [38]. These findings suggested that similar salt tolerance mechanism could exist in cowpea. Our candidate gene search based on trait-associated SNPs supports the involvement of vacuolar proteins in salt tolerance in cowpea. Vacuolar proteins are critical in maintaining the trans-Golgi network (TGN) during salt tolerance [39]. In addition, Vigun03g290600.1 has been reported to encode for xyloglucan:xyloglucosyl transferase that is induced upon salt stress in Arabidopsis thaliana [40]. Despite of the possible relationship existing between functional annotations of the candidate genes identified in this study and salt tolerance mechanism, further studies including transcriptomic analysis would be required to validate these genes.
Soil salinity has been shown to be a growing threat to agriculture worldwide [9]. Cowpea can be significantly impaired by soil salinity [8]. In this investigation we reported the variation of salt tolerance in a MAGIC cowpea population and identified salt-tolerant MAGIC lines. This MAGIC population has been registered [28], but lacking information on tolerance to salt stress. Therefore, our results complement the previously collected information on this MAGIC population and helps further increase its usefulness in cowpea breeding.

Conclusions
Our findings contribute to providing useful information for basic and applied research in genetic improvement to address soil salinity. Large variation in salt tolerance among the cowpea MAGIC lines has been identified. The salt-tolerant lines could be used as donor parents in breeding for salt tolerance in cowpea. In addition, a large number of cowpea SNP markers were found within or in the vicinity of genes that were reported to be directly involved in salt tolerance in other crops. Therefore, these SNPs can provide a framework for map-based cloning and designing gene-based markers in cowpea and other crops based on synteny. Breeding for salt tolerance via marker-assisted selection (MAS) would also be possible. However, additional studies would be needed to validate the candidate genes identified in this study.

Population development and genotyping
The MAGIC population development was previously described [28]. In brief, the first crosses between eight founder parents (IT89KD-288, IT84S-2049, CB27, IT82E-18, Suvita 2, IT00K-1263, IT84S-2246, and IT93K-503-1) were conducted in 2011. These parents were elite cultivars and breeding lines from Sub-Saharan Africa and the United States [41][42][43][44][45][46]. Greenhouse-grown seeds of 305 F 8:10 RIL lines and eight founders were obtained from University of California, Riverside (UCR). Ten seeds per line were hand-planted in a 5-foot-long row at the research station of the University of Arkansas, Fayetteville during the summer of 2018. Some lines were not able to flower due to photoperiodism under the Arkansas climate. A total of 234 MAGIC lines and eight parents were harvested. Seeds from each row were bulked for use in the salt tolerance evaluation.
The MAGIC population and the founders were genotyped using a total of 51,128 SNPs obtained from the Illumina Cowpea Consortium Array [47]. The SNP data as well as genetic diversity analysis of this population were previously reported [28]. After filtering, a total of 32,047 SNPs (missing data<10%, heterozygosity<10%, and minor allele frequency>5%) were used for further analysis.

Growth conditions and experimental design
Salt tolerance evaluation was conducted using a previously described methodology [48]. The experiment was carried out in the greenhouse of Harry R. Rosen Alternative Pest Control of the University of Arkansas, Fayetteville where the average temperature was 26 °C/21 °C (day/light) and the daylight length was 14 h (Fig. 7). Cowpea seeds were sown in pots previously filled up with 100 g Sunshine Natural & Organic (Agawam, MA). At the bottom of each pot was designed holes to prevent water from being stored during irrigation, which could lead to plant root asphyxia. In addition, paper towels were established at the bottom of each pot to prevent soil from leaking during irrigation. In each pot, a total of 8 seeds were sown and thinned to a total of 4 vigorous and uniform plants at one week after emergence. Plants were fertilized weekly by applying a solution of 50 mL of Miracle-Gro fertilizers (Scotts Miracle-Gro, Detroit, MI) to each pot.
The experiment was a randomized complete block design (RCBD) with two blocks and two replications within each block. Pots containing cowpea plants were placed on rectangular plastic trays to make the irrigation process convenient. For each line per replication, one pot was used as control receiving deionized water during irrigation, and the other pot receiving a salt treatment, so there were one salt-treated and one non-salt-treated pot for each replication.
The salt treatment (NaCl) started when the first trifoliate leaf began to expand (V1 stage) [49]. Salt concentration was 200 mM NaCl as previously suggested [23,[50][51][52]. Irrigation was applied to each tray containing 12 pots with either deionized water or salt solution. Irrigation was maintained in such a way that two-third of pot height was soaked with the treatment solution.
One salt-tolerant cowpea genotype ('09-529') and one salt-susceptible cowpea genotype (PI255774) were also included as checks [21,23]. Two subsets of the extremes, including ten most salt-tolerant and ten most salt-susceptible lines, were repeated at the end of the experiments.

In vivo chlorophyll measurement
Leaf chlorophyll was measured using a SPAD-502 Plus (Spectrum Technologies, Inc., Plainfield, IL). Measurements were achieved at one day prior to salt treatment and when the susceptible controls were completely dead, which was about 14 days after the first salt stress. For each plant, chlorophyll measurement was conducted three times on both trifoliate and unifoliate times, respectively, and the average read was recorded and analyzed. Measurements were done on three different positions on the leaf surface in order to limit the edge effect [21,23]. Data were taken from all plants under salt stress and non-salt stress conditions.

Plant fresh above-ground parameter
Plant height was recorded for each cowpea seedling at one day before the salt treatment began and when the susceptible controls were dead, indicative of the end of plant growth in the susceptible genotypes [23]. Data on both fresh leaf and fresh stem biomass from each plant were also taken. The above-ground fresh biomass corresponded to the sum of fresh leaf and stem weights.

Leaf injury scoring
Leaf injury score has been successfully used as a reliable parameter for screening salt tolerance at the seedling stage in cowpea [23]. It has been shown to be highly correlated with Na + and Cl − contents in leaves [53], and can accurately assess salt tolerance/susceptibility when leaf ion extraction is costly [23,53]. Leaf injury scoring was conducted when the susceptible controls were completely dead, using a previously established scale (1 = healthy plants, 2 = sign of leaf chlorosis, 3 = expansion of chlorosis on leaf surface, 4 = totally chlorotic leaf, 5 = first sign of necrosis, 6 = expansion of necrosis on leaf surface, and 7 = completely dead plants) .

Phenotypic data analysis
Relative tolerance indexes (RTIs) for chlorophyll, plant height, fresh leaf biomass, fresh stem biomass, and total fresh above-ground biomass were used to assess the impact of salt stress relative to the non-salt stress condition. RTIs were calculated as following [23,54]. Data distribution was visualized using the MASS package of the R ® 3.6.1. version. The analysis of variance (ANOVA) was done using PROC MIXED of SAS ® 9.4 (SAS Institute Inc., Cary, NC). Mean separation was conducted using a protected least significant difference (LSD) procedure at α = 0.05. LSD procedure was defined as LSD=t α/2 √2MSError/n, with t α/2 being the critical value from the t-table and having a degree of freedom [df(SSError)] corresponding to the difference between the number of observations and the number of replications, and n being the number of replications. The ANOVA model was as follows.
with µ being the overall mean, Y i(j)k being the response from the k th genotype (G k ) (fixed effect) at the i th replication (R i(j) ), which was nested under the j th run (block) (T j )(fixed effect), and TG jk being the interaction effect between the k th genotype (G k ) and the j th run (block) (T j ).
The broad sense heritability (H) was estimated using the following formula [55].
with σ 2 G being the total genetic variance, σ 2 GXR being the Genotype X Run variance, σ 2 e being the residual variance, n b being the number of runs, and n r being the number of replications. The estimates for σ 2 G and σ 2 GXR were [EMS(G)-EMS(GXB)]/ n b *n r and [EMS(GXB)-Var(Residual)]/n r . EMS(G), EMS(GXB), and Var(Residual) were obtained from the ANOVA table.
Pearson's correlation coefficients between the average number of dead plants per pot, average leaf injury score, fresh leaf biomass under salt stress, SPAD chlorophyll content under salt stress, RTI_C, RTI_H, RTI_FL, RTI_FS, and RTI_FB were calculated using the R ® 3.6.1. version. A chord diagram was used in order to better visualize the pairwise correlation between traits. Chord diagram was established in R ® v.3.6.1 using the package 'circlize' [56].

Genome-wide association study (GWAS)
GWAS was conducted using a Bayesian Information and Linkage Disequilibrium Iteratively Nested Keyway (BLINK) model and run in R ® 3.6.1 using the package 'BLINK' [57]. BLINK has been demonstrated to have an enhanced statistical power and to be more efficient compared to previously developed models [57]. LOD threshold was set to 3 [58].
The two FEM models in BLINK were defined as following. Where y: being the vector phenotype; M i1 , M i2 , …, M ik : the genotypes of k pseudo QTNs that were initially empty and with effects b 1 , b 2 , …, b k , respectively; M ij : the j th genetic marker of the i th sample; and e i : the residual Y i(j)k = µ+T j + G k + R i(j) + TG jk +ε i(j)k where i = 1, 2, j = 1, 2, and k = 1 . . . 234 having a distribution with mean zero and a variance σ 2 e . LD heatmaps were generated using the package 'LDheatmap' in R ® 3.6.1 [59]. Overlapping SNP markers between different traits were visualized using a Venn diagram that was established using the online software program accessible at http:// jvenn. toulo use. inra. fr/ app/ examp le. html.

Candidate gene discovery
A 40-kb genomic region harboring each significant SNP was used for searching candidate genes in the Phytozome 12 database (https:// phyto zome. jgi. doe. gov/), with an emphasis on those with functional annotations relevant to abiotic stresses.

Epistatic interaction modelling
Pairwise epistatic interaction analysis (SNP * SNP interaction) was conducted using PLINK v1.07 [60]. The command line for conducting epistasis analysis in PLINK was 'plink --file mydata --epistasis' . The interaction effect of two SNPs was estimated using the following model [60].
with E[Y|Snp i , Snp j ] being the vector of expected values for the response given the SNP data, β 0 being the intercept, β i being the main effect for the Snp i , β j being the main effect for the Snp j , and β ij being the interaction effect (epistasis) between Snp i and Snp j . The parameter of interest in the above model was β ij and the test to be conducted was H0: β ij = 0. Choosing a minimum p-value for declaring a significant interaction effect can inflate the Type 1 error rate [61]. However, the current approach using various techniques for identifying a significant threshold while reducing the bias in estimating β ij and limiting the Type 1 error rate could be still extremely computationally intensive. Therefore, we used an arbitrary threshold (p-value ≤ 10 − 6 ) in this study given the great number of possible pairwise interactions and practical reasons during the data visualization process, while being biologically reasonable. Pairwise epistatic interaction was visualized using the package 'circlize' and run in R ® 3.6.1 [56].

Statement on Research involving plants
The study protocol complies with relevant institutional, national, and international guidelines and legislation.