GWAS and WGCNA uncover hub genes controlling salt tolerance in maize (Zea mays L.) seedlings

Two hub genesGRMZM2G075104andGRMZM2G333183involved in salt tolerance were identified by GWAS and WGCNA. Furthermore, they were verified to affect salt tolerance by candidate gene association analysis. Salt stress influences maize growth and development. To decode the genetic basis and hub genes controlling salt tolerance is a meaningful exploration for cultivating salt-tolerant maize varieties. Herein, we used an association panel consisting of 305 lines to identify the genetic loci responsible for Na+- and K+-related traits in maize seedlings. Under the salt stress, seven significant single nucleotide polymorphisms were identified using a genome-wide association study, and 120 genes were obtained by scanning the linkage disequilibrium regions of these loci. According to the transcriptome data of the above 120 genes under salinity treatment, we conducted a weighted gene co-expression network analysis. Combined the gene annotations, two SNaC/SKC (shoot Na+ content/shoot K+ content)-associated genes GRMZM2G075104 and GRMZM2G333183 were finally identified as the hub genes involved in salt tolerance. Subsequently, these two genes were verified to affect salt tolerance of maize seedlings by candidate gene association analysis. Haplotypes TTGTCCG-CT and CTT were determined as favorable/salt-tolerance haplotypes for GRMZM2G075104 and GRMZM2G333183, respectively. These findings provide novel insights into genetic architectures underlying maize salt tolerance and contribute to the cultivation of salt-tolerant varieties in maize.


Introduction
Salinity is a major abiotic stress affecting crop growth, development, and yield (Cui et al. 2015;Zhu 2001). Globally, approximately 50% of irrigated lands and 20% cultivated lands are affected by salt stress (Hu and Schmidhalter 2005;Luo et al. 2019a). As an important crop, maize is sensitive to salt stress, particularly at the seedling stage (Cui et al. 2015). Understanding the genetic architecture and causal genes of salt tolerance in maize seedlings is an urgent need to reduce salinity damage and cultivate salt-tolerant maize cultivars for breeders.
Two important mechanisms underlying the response of salinity stress have been widely verified in plants, namely SOS (salt overly sensitive) and HKT (high-affinity K + transporter) pathways (Hamamoto et al. 2015;Zhu 2002). In SOS pathway, three genes SOS1, SOS2, and SOS3 jointly participate in the plant metabolism under salt stress (Zhu 2002). SOS1 is a Na + /K + antiporter, which specifically expressed in root tips, and epidermal cells of stems and leaves (Shi et al. 2000(Shi et al. , 2002. Under salinity stress, SOS1 transports Na + out of cells via H + electrochemical potential (Shi et al. 2002). SOS2 and SOS3 encode Serine/Threonine protein kinase and calcium-binding protein with EF hand motifs, respectively Liu et al. 2000). During salt stress, Ca 2+ activates SOS2 and causes the formation of SOS2-SOS3 complex, which then activates the C-terminal region of SOS1 by phosphorylation (Quintero et al. 2011). Proteins MPK6 and SCaBP also participate in the SOS pathway (Ji Communicated by Benjamin Stich. 1 3 et al. 2013). HKT is an Na + -preferential transporter, which mainly controls root-shoot Na + delivery after the Na + transits xylem and gets into the cell (Munns 2002;Ren et al. 2005;Zhang et al. 2018). Additionally, IAA (auxin), ABA (abscisic acid), ETH (ethylene) and other plant hormones were reported to be the key regulators in the responses of salt stress (Ryu and Cho 2015).
Linkage mapping is a useful method to dissect the genetic basis of quantitative traits (Cui et al. 2015). Currently, numerous quantitative trait loci (QTL) controlling salt tolerance have been identified in crops by linkage mapping. In rice, SKC1 was the first cloned QTL for salt tolerance in rice (Ren et al. 2005). Using a RIL (F 2:9 ) population detected 16 QTL associated with rice salt-tolerance (Wang et al. 2011). Hamwieh and Xu (2008) uncovered a major salt-tolerant QTL, which explained 68.7% of phenotypic variation in Soybean. In maize, Hoque (2013) identified 10 salt tolerance QTL in an F 2 population consisting of 163 lines. QTL analyses using an F 2:5 population revealed eight QTL related to salt tolerance-related traits (Cui et al. 2015). Recently, Zhang et al. (2018) cloned a salt tolerance-associated gene ZmHKT1 based on QTL fine-mapping in maize. However, only a few QTL (SKC1 in rice and ZmHKT1 in maize) were fine-mapped and cloned in crops, owing to the large confidence intervals of these QTL (Ren et al. 2005;Zhang et al. 2018). Genome-wide association studies (GWAS) provide the detection of QTL with a high-resolution approach (Sukumaran et al., 2015) and have been extensively used for excavating the candidate genes controlling targeted traits (Korte and Farlow 2013). Using GWAS, 66 significant single nucleotide polymorphisms (SNPs) were found to be associated with 12 salt tolerance traits in rice (Kumar et al. 2015). In barley, GWAS using 2,671 lines detected a salt tolerance-associated gene HKT1;5 (Hazzouri et al. 2018). In maize, a Na + transport-associated gene ZmHAK4 was identified via GWAS and subsequently functionally-validated by overexpression and mutant strategies . Recently, Luo et al. (2021) detected two genes ZmCLCg and ZmPMP3 using GWAS, which were confirmed to regulate the maize salt tolerance (Luo et al. 2021). Weighted gene co-expression network analysis (WGCNA) facilitates the revelation of the core gene networks based on the gene expression patterns from RNA-seq data, which has become a favored and popular technique in discovering hub factors controlling traits (Wan et al. 2018;Guo et al. 2020). Liu et al. (2016a) identified six priority transcription factors associated with lysine biosynthesis in the LBPGs module of WGCNA results. Two hub genes involved in element accumulations were detected using WGCNA and further verified by a mutant-based approach (Schaefer et al. 2018).
In plants, salt tolerance was previously evaluated by several traits, such as shoot fresh weight, shoot dry weight, plant height (Cui et al. 2015;Kumar et al. 2015). By contrast, Na + and K + accumulations are more direct and important traits for salt tolerance due to their important roles in regulating osmotic tolerance and Na + /K + delivery under salt stress ). In the present study, GWAS was performed in a panel consisting of 305 lines to detect significant SNPs correlating with Na + and K + accumulations under salinity conditions. Subsequently, expression values of these candidate genes were utilized to carry out WGCNA for identifying the priority genes. Finally, the target genes were individually PCR-amplified among 77 randomly selected lines from the association panel and the variation loci for each gene were then subjected to candidate gene association analysis. The goal of this study was to determine the hub candidate genes affecting Na + -and K + -related traits and identify the favorable haplotypes of the hub genes for improving the salt tolerance in maize breeding.

Plant materials and trait measurements
An association panel including 305 lines was used to carry out GWAS in this study. This panel was composed of tropical, stiff stalk (SS), non-stiff stalk (NSS), and other germplasms (Zhang et al. 2016).
A randomized complete block design was used for the trials with three replications for each line. Thirty kernels of each line were treated with 10% H 2 O 2 for 10 min and then soaked with saturated CaSO 4 for 12 h. After being rinsed three times by deionized water, the seeds were placed on vermiculites for germination in the growth chamber with a 16/8 h light/darkness cycle, a 25/22℃ day/night temperature cycle, a 200 μmol photonsm −2 s −1 brightness, and a relative humidity of 65% . At the three-leaf stage, the seedlings were divided into two groups and separately cultivated in Hoagland's solutions (normal treatment, NT) and Hoagland's solutions containing 150 mM NaCl (salt stress, SS) for 7d. The formulas of Hoagland's solution were described by Cui et al. (2015). In each line, eight seedlings with healthy and consistent growth for each treatment were selected and then kept in an oven at 80 °C for 72 h.
The dried shoots and roots were, respectively, crushed and approximately 2 g powder of per sample was diluted with 80% H 2 SO 4 for extracting Na + and K + by MARS6 instrument. Inductively Coupled Plasma Mass Spectrometer (ICP-MS) was used to measure Na + and K + contents.
The mean value, standard deviation (SD), maximum value, and minimum value of each trait were calculated by Microsoft Excel. SPSS software (Statistical Product and Service Solutions, version 21.0, IBM, Armonk, NY) was used to analyze the correlations of the collected traits and compare the phenotypes between two treatments (normal and salt treatments) for each trait by t test. Analysis of variance (ANOVA) was based on the model Here y ij , μ, R i , G j , and E ij denote the observation from the ijth plot, the overall mean, the effect of the ith replication, the effect of the jth genotype, and the experimental error, respectively ). SAS 9.3 (Statistical Analysis System, version 9.3, SAS Institute, Cary, NC, USA) was utilized to estimate the broad-sense heritability (H 2 ), which was described as follows (Pace et al. 2015): where, 2 G , 2 P , MSG, MSE, and rep represent genotype variance, phenotype variance, mean square of genotype, mean square of error, and the replication numbers (= 3), respectively.

Genotypic data analysis
Genotypic data of the 305 lines were obtained using the Illumina MaizeSNP50 BeadChip, as described in our previous study (Zhang et al. 2016). A total of 46, 408 SNPs evenly distributed over 10 chromosomes were remained after removing the SNPs with missing rate > 20%, heterozygosity > 20%, or minor allele frequency (MAF) < 0.05. In previous studies, evaluation of population structure and linkage disequilibrium (LD) in this panel revealed that the 305 lines can be divided into three subgroups and the LD decay approximately reached to 300 Kb (r 2 = 0.1) (Zhao et al. 2018;Zhang et al. 2020a). Thus, PCA (principal component analysis) = 3 was used as a covariate in GWAS, and the gene models located within the 300 Kb upstream and downstream of the significant SNPs were considered as the candidate genes in this study.

GWAS
Currently, three models were adopted in GWAS, namely GLM (General Linear Model), FarmCPU (Fixed and random model Circulating Probability Unification) model, and MLM (Mixed Linear Model). GLM trends to generate false positives due to no covariates included in this model, whereas MLM results in many false negatives due to the stringent redresses of population structure (Q matrix) and kinship (K matrix) (Hu et al. 2017;Ma et al. 2020). Farm-CPU model addresses the confounding problem of testing markers using iterative algorithms, which is moderately stringent in comparison to GLM and MLM (Liu et al. 2016b;Ma et al. 2020). Consequently, FarmCPU model was used in this study for GWAS. As a covariate, PCA was calculated by GAPIT software (Lipka et al. 2012) and then added to FarmCPU model. GAPIT software and FarmCPU model were both run in R studio (version 3.4.1). Herein, the significance threshold was set as: P-value = 0.05/n (n is the total SNP number) = 0.05/46,408 = 1.08 × 10 -6 . SNPs with P-value < 1.08 × 10 -6 were considered as the significant SNP-trait associations. Moreover, the genes located in the LD regions of the significant SNPs were identified as candidate genes affecting traits.

WGCNA and hub gene identification
The gene expression values were obtained from RNA-seq in our previous study (Accession number: CRA003872; Zhang et al. 2021). Briefly, the seeds of L2010-3 (a salt-tolerant line) and BML1234 (a salt-sensitive line) were germinated in a growth chamber. At two-leaf stage, the seedlings were divided into two groups and, respectively, transplanted into Hoagland's solutions (control) and Hoagland's solutions with 150 mM NaCl (salt treatment) to culture in a growth room. The roots from L2010-3 and BML1234 were separately collected at 0 h (control), 6 h (control, treatment), 18 h (control, treatment), and 36 h (control, treatment) with two biological replicates, and a total of 28 samples were subjected to transcriptome sequencing. The obtained raw sequencing reads were filtered by fastp software and the clean reads were retained. Gene expression levels were calculated by the clean reads mapped on B73 genome and were then normalized to transcripts per kilobase million (TPM) via consumer script. The TPM of candidate genes detected by GWAS were subjected to WGCNA based on the WGCNA package in R studio software (Langfelder and Horvath 2008). The parameters of WGCNA program were as follows: variance data expression > 0; no missing data expression < 0.1; soft threshold = 8 (estimate value); max block size = 200; deep split = 2; min module size = 3; merge cut height = 0.1. In each module, the genes with eigengenebased connectivity value (|KME|) > 0.9 and topological overlap measure (TOM) value > 0.2 were regarded as hub genes (Buckberry et al. 2017;Zhang and Horvath 2005). Function annotations of the candidate genes were obtained by searching NCBI website (RefGen_v2) (https:// www. ncbi. nlm. nih. gov/). Kyoto encyclopedia of genes and genomes (KEGG) pathway analyses were implemented in GENE DENOVO platform (https:// www. omics hare. com/ tools/ Home/ Soft/ pathw aygsea). Cytoscape software was used to draw the regulation networks of target genes (Shannon et al. 2003).

Candidate gene association study
Among these hub genes, the ones that were reported to correlate with salt tolerance of plant species were considered as prioritized hub genes and were individually subjected to candidate gene association study. The gene body and its upstream 2,000 bp sequences were PCR-amplified in 77 randomly selected inbred lines from the association panel. Primers were designed in Primer 3.0 (v. 0.4.0). DNAMAN software (version5.2.2, lynnon bio-soft, Canada) was used to perform sequence alignment between the amplified sequences and B73 genome sequence to identify sequence diversity including SNPs and insertion/deletion (InDels) (Zhang et al. 2020b). The associations between the phenotypes and SNPs/InDels with MAF ≥ 0.05 were analyzed by conducting candidate gene association based on a GLM + PCA model in TASSEL 4.0 software with the significant threshold P = 0.01 (Bradbury et al. 2007). LD decay between the SNPs was calculated by HaploView software (http:// www. broad. mit. edu/ mpg/ haplo view). Haplotypes were identified by using significant markers located in target genes.

Superior allele analysis in elite lines
For each significant SNP detected in SS conditions, a superior allele was determined by the effect value. For each SNP, SP (%) = (NS/TN) × 100%, SP is the superior allele proportion, NS and TN represent the number of superior alleles in the elite lines and the total number of alleles in the elite lines, respectively (Ma et al. 2018). Additionally, the heatmap of superior allele distributions was drawn in the R software using a heatmap package.

Evaluation of salt-tolerance phenotypes
To systematically evaluate the phenotypes of salt-tolerance, mean value, standard deviation (SD), maximum value, and minimum value were collected for each trait. For all traits, highly phenotypic variations were observed with SDs ranged from 0.12 to 937.81 (Table S1). Significant differences (P < 0.001) for all traits were detected between NT and SS ( Fig. 1), suggesting salt stress in the association panel was effective. Under SS conditions, SNaC, RNaC, SNaC/SKC, RNaC/RKC, and KTC were obviously enhanced with the increased folds ranged from 0.54 to 18.32, whereas SKC, RKC, and NaTC were greatly lowered with the decreased folds ranged from 0.26 to 0.52 in comparison to their corresponding phenotypes under NT condition ( Fig. 1; Table S1). Under NT and SS conditions, the phenotype frequency distributions of all traits displayed normal distributions (Fig. 1). Broad-sense heritability (H B 2 ) estimates of traits were in a range of 0.35-0.89 across two conditions (Table S1). Under NT condition, SNaC and SNaC/SKC displayed the highest positive correlation (r = 0.88), whereas RKC and KTC showed the highest negative correlation (r= −0.65) (Table S2). Under SS condition, the highest positive correlation was observed between NaTC and SNaC, with r = 0.80, and the highest negative correlation was observed between SNaC and SNaC/SKC, with r = −0.63 (Table S2).

WGCNA
To further identify the hub genes out of all 120 candidate genes from GWAS, we conducted a WGCNA using the expression values (TPM) of these genes. As a result, these 120 genes were classed into four subgroups, namely blue, brown, turquoise, and yellow modules (Fig. 4a). Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis revealed that the genes Fig. 3 Significant SNPs detected in SS condition. SS represents salt stress. a Manhattan plot of significant SNPs. The numerals i, ii, and iii represent KTC (K + transport coefficient), RKC (root K + content), and SNaC/SKC (shoot Na + content/shoot K + content), respectively.
Red dots represent significant SNPs associated with target traits. b, c, d QQ plots of significant SNPs. SNaC/SKC, RKC, KTC represent shoot Na + content/shoot K + content, root K + content, and K + transport coefficient, respectively 1 3 in the blue, brown, and turquoise modules were significantly (P < 0.05) enriched in five (glucosinolate biosynthesis, valine, leucine and isoleucine biosynthesis, pantothenate and CoA biosynthesis, ribosome, and valine, leucine and isoleucine degradation), two (RNA transport, and protein processing in endoplasmic reticulum), one (ABC transporters) pathways, respectively. As for the yellow module, no significant pathway was enriched. Additionally, a total of nine genes were discovered and subsequently defined as hub genes with the principles of |KME|> 0.9 and TOM > 0.2 (Table 2) (Table 2). Most strikingly, ABC transporters were previously reported to be key factors in regulating ion exchange under salt stress, therefore, turquoise module was focused in our study (Moons, 2003;Saha et al. 2015). In turquoise module, GRMZM2G075104 and GRMZM2G333183, respectively, annotated as ubiquitin conjugation factor E4 and ABC transporter B family member 14 (Table 2), which participates in salt stress response in  . As such, GRMZM2G075104 and GRMZM2G333183 were considered as prioritized hub genes for further study. According to the node cytoscape results from WGCNA, we constructed the regulation networks of the genes located in turquoise module (Fig. 4b).

Distribution of superior alleles in elite inbred lines
Since our association panel included 30 elite lines which were collected from maize breeding programs (Ma et al. 2018), this allowed us to estimate the utilization of superior alleles during maize breeding. According to the GWAS results, significant SNPs were detected for SNaC/SKC, RKC, and KTC under salt stress condition. Herein, the allele associated with a lower phenotype was defined as the superior allele for SNaC/SKC, whereas the allele with a higher phenotype was designated as the superior allele for RKC and KTC. Among the 30 elite inbred lines, the superior allele percentages of significant SNPs ranged from 10.71% to 93.00%, with five of SNPs (SYN38412, PZE-104132392, SYN6394, PZE-110042260, and SYN24195) containing ≥ 50% superior alleles, whereas the remaining two SNPs (SYN26336 and PUT-163a-21390620-1531) containing < 50% superior alleles (Fig. 7). Moreover, 20 elite lines that contained 4-7 superior alleles displayed higher phenotypic values, with increased percentages of 22.92% (RKC), and 65.91% (KTC), relative to the other 10 elite lines that contained 0-3 superior alleles. However, for SNaC/SKC, the above 20 lines that contained 4-7 superior alleles had a lower mean value in comparison to the other 10 lines that contained 0-3 superior alleles, with the reduced percentage of 55.52%. These findings suggested that the superior alleles have obviously additive effects on salt tolerance.

Using Na + and K + -related traits for evaluating salt tolerance in maize
In previous studies, seedling shoot length, root length, shoot fresh weight, root fresh weight, shoot dry weight, root dry weight, tissue water content, plant height, and plant survival rate were used to evaluate the salt tolerance in maize (Cui et al. 2015;Luo et al. 2017Luo et al. , 2019aLuo et al. , b, 2021. According to the genetic architecture analysis of these traits, several QTL controlling maize salt tolerance were identified. However, all these traits represent indirect responses to salt stress. Under salt stress, Na + and K + -related traits (SNaC, RNaC, SKC, RKC, SNaC/ SKC, RNaC/RKC, KTC, and NaTC) are direct evaluation standards for salt tolerance of plants. Meanwhile, these traits are involved in two salt-plant signals, SOS and HKT pathways. Recently, ZmHKT1 (Na + transporter), ZmHAK4 (Na + transporter) and ZmHKT2 (K + transporter) were demonstrated to mediate salt tolerance of maize under salinity conditions (Cao et al. 2019;Zhang et al. 2018Zhang et al. , 2019. Therefore, Na + and K + -related traits were taken as evaluation indexes for salt tolerance of maize seedlings in this study. Abundant phenotypic variations were observed for all measured traits in this association panel, with the variation coefficients of 24.01-50.00% and 27.21-89.24% in NT and SS, respectively (Table S1), suggesting this panel was suitable for GWAS of the eight Fig. 5 Candidate gene association analysis of GRMZM2G075104. a Significant markers associated with SKC (shoot K + content) and SNaC/SKC (shoot Na + content/shoot K + content). b Pairwise LDs (linkage disequilibriums) between the markers. c Comparison of SKC (shoot K + content) between HapI (haplotype I) and HapII (haplotype II). *Significant at P < 0.05. d Comparison of SNaC/SKC (shoot Na + content/shoot K + content) between HapI (haplotype I) and HapII (haplotype II). **Significant at P < 0.01. e Details of HapI (haplotype I) and HapII (haplotype II). "-" Represents insertion or deletion. N represents the inbred line numbers of each haplotype Na + and K + -related traits under two conditions. Herein, the heritability values of the traits ranged from 0.35 to 0.89. Under SS condition, seven SNPs were detected to control SNaC/SKC, RKC, and KTC, with the total PVE ranging from 6.14 to 72.30% among these traits. The unexplained PVE probably resulted from (i) the minoreffect SNPs which were missing; (ii) the lower SNP density of the panel used in this study.

A combination of GWAS and WGCNA is an efficient strategy for elucidating genetic architecture and hub genes of salt tolerance in maize
Currently, GWAS is widely used in dissecting genetic basis of complex agronomic traits. Using GWAS, previous studies revealed many candidate genes controlling target traits in maize. For embryonic callus regenerative capacity in maize, a total of 40 candidate genes were identified using GWAS, including a previously verified gene WOX2 (Ma et al. 2018). Zhang et al. (2020a) found an ear row number-associated SNP SYN8062 in an association panel, which was closely linked to the gene model GRMZM2G098557. Subsequently, GRMZM2G098557 was validated to regulate ear row number by using a segregation population. Similarly, several genes associated with salt tolerance were also detected by GWAS in maize, such as ZmCLCg, ZmPMP3, and ZmHAK4 (Luo et al. 2021;Zhang et al. 2019). In this study, we totally identified 120 candidate genes using GWAS, 65 of which were located in the previously detected QTL for maize salt tolerance. Specifically, a total of 20, 25, 2, and 18 genes were situated in two (qFL1 and qSL1), one (qNPH9-2), one (qPHI3), and five (qNPH9-2, qRF9, qSF9, qFF9, and qRD9) QTL, respectively (Luo et al. 2017(Luo et al. , 2019a. Among the candidate genes, GRMZM2G104783 (PZE-104132392), involved in SNaC/SKC was annotated as a cytochrome P450 superfamily protein (Table S4). In Arabidopsis thaliana, a cyp709b3 (a member of cytochrome P450 superfamily) mutant showed sensitivity to salt during germination and low tolerance to salt stress (NaCl concentration: 150 mM) (Mao et al. 2013). In Gossypium hirsutum, two cytochrome P450 members, Gh_D07G1197 and Gh_A13G2057 played positive roles in enhancing drought and salt stress tolerance (Magwanga et al. 2019). A SNaC/ SKC-associated gene, GRMZM2G372952 (SYN6394) was annotated as a zinc finger protein ZAT4 (Table S4). Sun et al. (2019) reported that the salt and drought tolerances were enhanced by overexpressing the GMZAT4 in Arabidopsis thaliana. The SNaC/SKC-correlated gene GRMZM2G372632 (SYN6394) encodes a pentatricopeptide repeat-containing protein At4g31850 in chloroplastic (Table S4), the homologous gene of which, wsl was reported to enhance salinity sensitivity in rice (Tan et al. 2014). GRMZM2G171428 (SYN6394) that was associated with SNaC/SKC encodes a WRKY transcription factor 4 (Table S4). Previous studies showed that a WRKY transcription factor WRKY17 negatively regulated salt stress tolerance in transgenic Nicotiana benthamiana and Arabidopsis via abscisic acid signaling pathway (Yan et al. 2014;Cai et al. 2017). However, overexpression of another WRKY transcription factor VvWRKY30 in Arabidopsis increased its resistance to salt stress at different growth stages (Zhu et al. 2019). A RKC-associated gene GRMZM2G378852 (SYN26336) encodes a MAPKKK family protein kinase (Table S4), whose homologous gene, At1g73660 was demonstrated to negatively regulate the salt tolerance in Arabidopsis (Gao and Xiang 2008). Gene model GRMZM2G070304 (PUT-163a-21390620-1531) that was involved in KTC was annotated as glycerol-3-phosphate acyltransferase 3 (GPAT3) ( Table S4). Sui et al. (2017) confirmed that the salt tolerance of Arabidopsis could be improved by overexpressing the homologous gene (SsGPAT) of GRMZM2G070304 from Suaeda salsa. Moreover, two genes GRMZM2G075104 and GRMZM2G333183, belonging to ubiquitin conjugation factor E4 and ABC transporter B family member 14, respectively (Table S4), were reported to be associated with salt stress tolerance in Spinach (Spinacia oleracea L.), Arabidopsis, Oryza sativa and other species (Bagheri et al. 2015;Lee and Kim 2011;Moons 2003;Saha et al. 2015).
Recently, to decode the genetic control and hub candidate genes related to target traits, a joint analysis of GWAS and WGCNA has become an efficient approach. Using this strategy, two hub genes GRMZM2G477685 and GRMZM2G135536 controlling primary root growth were identified by Hwang et al. (2018). Similarly, Guo et al. (2020) detected seven hub candidate genes involved in seminal root length of maize seedlings under drought stress though GWAS and WGCNA. In our study, we totally identified nine candidate genes associated with salt tolerance in maize seedlings by using a combination of GWAS and WGCNA (Table 2). Combined with the gene annotations, GRMZM2G075104 and GRMZM2G333183 were finally determined as hub genes. Therefore, these two promising genes were prioritized for further study.

Interaction predictions of hub genes
In this study, the hub genes GRMZM2G075104 and GRMZM2G333183 were both located in turquoise module. To identify the upstream regulators of the two hub genes, we further focused on these transcript factors in turquoise module. Among the 19 genes co-expressed with GRMZM2G333183 (Fig. 4b), only GRMZM2G171428 was annotated as transcript factor (WRKY transcription factor 4) in this module. Interestingly, we found a WRKY binding motif (TAG TCA AA) in the promoter region (upstream 3,000 bp) of GRMZM2G333183 (P value = 9.62E-05) using the plant transcriptional regulatory map website (http:// plant regmap. gao-lab. org/ index. php). This finding suggested that the WRKY transcription factor, GRMZM2G171428 probably regulate the hub gene GRMZM2G333183 by binding the motif (TAG TCA AA) in the promoter of GRMZM2G333183. Utilizing STRING software (https:// string-db. org/), we subsequently predicted the encoding protein interactions between the hub genes and the other genes in the co-expression network. Eventually, we did not search out any interaction in turquoise module for hub genes GRMZM2G075104 and GRMZM2G333183, indicating the encoding proteins of hub genes performed functions independently.

Identification of superior haplotypes for hub genes
Currently, traditional breeding methods in crops produce lower genetic gains than required to achieve food sufficiency. In this context, advanced breeding technologies should be proposed and utilized. Superior haplotypes affecting many agronomic traits (such as yield, kernel quality, and abiotic stress responsiveness) have been verified across maize (Yu et al. 2019;Li et al. 2020;Zhang et al. 2020b). Genetic variations occurred in promoter and coding regions probably leads to the variations of gene expression abundances and protein sequences. Candidate gene association analysis was able to identify the significant variations based on a P-value threshold and thus widely used in detecting/dividing superior haplotypes of target genes. Using this method, Li et al. (2020) analyzed the haplotypes of candidate genes associated with trait ear tip-barrenness in maize. Zhang et al. (2020b) reported one favorable haplotype for gene Zm00001d047868 which influenced maize kernel moisture content. Additionally, a previous study reported that overexpressing favorable haplotype of ZmEREB180 increased maize waterlogging tolerance (Yu et al. 2019). In this study, the favorable haplotypes of GRMZM2G075104 and GRMZM2G333183 were determined as TTG TCC G-CT and CTT using candidate gene association analysis, respectively. In future studies, to verify the functions of GRMZM2G075104 and GRMZM2G333183, the favorable haplotypes TTG TCC G-CT and CTT should be given priority for overexpression/knockout. In addition, superior haplotypes TTG TCC G-CT and CTT identified in this study can be used for molecular marker-assisted breeding of increasing salt tolerance in maize.

Application of superior alleles in cultivating salt-tolerance maize varieties
According to the distributions of seven significant SNPs uncovered in 30 elite inbred lines, we found two SNPs (SYN26336 and PUT-163a-21390620-1531) containing < 50% superior alleles (Fig. 7), which suggested that these alleles were not effectively selected during artificial selections. The possible reason is that maize salt tolerance was not a main trait concerned by breeders. The other five SNPs (SYN38412, PZE-104132392, SYN6394, PZE-110042260, and SYN24195) contained ≥ 50% superior alleles across the 30 elite lines (Fig. 7). Especially for PZE-104132392 and PZE-110042260, 27 superior alleles (90%) were observed among the 30 elite lines (Fig. 7). These findings implied that these alleles might linked with some major traits and were simultaneously maintained with the concerned traits in breeding programs. In further studies, the SNPs with low proportions of superior alleles should be emphasized by marker-assisted selection in cultivating salt-tolerance maize varieties. In addition, the inbred lines K22 and F06 contained seven (100%) and six (85.7%) superior alleles, respectively (Fig. 7), indicating these were salttolerance lines. The present findings implied that K22 and F06 can be used for gene functional research of salt stress and breeding of salt-tolerance maize varieties.