Genetic Variation in ZmPAT7 Contributes to Tassel Branch Number in Maize

Tassel branch number (TBN) is one of the important agronomic traits that contribute to the efficiency of seed production and has been selected strongly during the modern maize breeding process. However, the genetic mechanisms of TBN in maize are not entirely clear. In this study, we used a B73 × CML247 recombination inbred lines (RILs) population to detect quantitative trait loci (QTLs) for TBN. A total of four QTLs (qTBN2a, qTBN2b, qTBN4, and qTBN6) and six candidate genes were identified through expression analysis. Further, one of the candidates (GRMZM2G010011, ZmPAT7) encoding an S-acyltransferase was selected to validate its function by CRISPR-Cas9 technology, and its loss-of-function lines showed a significant increase in TBN. A key SNP(−101) variation in the promoter of ZmPAT7 was significantly associated with TBN. A total of 17 distant eQTLs associated with the expression of ZmPAT7 were identified in expression quantitative trait loci (eQTL) analysis, and ZmNAC3 may be a major factor involved in regulating ZmPAT7. These findings of the present study promote our understanding of the genetic basis of tassel architecture and provide new gene resources for maize breeding improvement.


Introduction
Maize is a monoecious cross-pollination crop. The male inflorescence (tassel) provides pollen for hybridization, which is important for maize grain production. Tassel size and weight have been reduced during modern breeding improvement [1,2]. Smaller tassels reduce the plant's energy requirements for general maintenance and make more energy available for grain production [3]. On the other hand, tassel with sufficient pollens directly contributes to the efficiency of seed production and hybrid yield performance [4]. Therefore, an optimal tassel architecture is needed that can produce sufficient pollen and transform more energy into the kernel. As a major component of tassel architecture, tassel branch number (TBN) is significantly negatively correlated with grain yield, and the reduction in TBN has been experienced strong selection during the modern maize breeding process [5].
TBN is a complex quantitative trait controlled by multigenes. A large number of quantitative trait loci (QTLs)/quantitative trait nucleotides (QTNs) for TBN have been mapped in different mapping populations. Zhao et al. (2017) detected 11 QTLs using two F 2:3 populations and identified 13 meta-QTLs (mQTLs) by integrating QTLs reported previously [6]. Li et al. (2018) summarized 100 QTLs and 360 QTNs for TBN, identified in several previous studies, and found that these QTLs and QTNs were mainly located at 20 genomic regions [7]. Recently, Wang et al. (2018) detected six minor effect QTLs for TBN using a Zheng58 × Chang7-2 recombination inbred lines (RILs) population with super-high-density genotyping [8]. Liu et al. (2019) detected seven QTLs using an F 2:3 population that generated from the cross between Lv28 and H082183 [9]. Wang et al. (2019) detected 51 QTNs and 19 QTLs associated with TBN in an association mapping population and a doubled haploid (DH) population, respectively [10]. Notably, a nested association mapping (NAM) population was used to identify 39 QTLs and 325 QTNs for TBN by joint linkage analysis and genome-wide association studies (GWAS) [11]. Furthermore, Wu et al. (2016) integrated two NAM populations (a NAM from the US and a NAM from China) and identified 63 QTLs and 549 QTNs for TBN [12]. Identification of those QTLs is beneficial to elucidating the genetic basis and facilitating gene cloning for TBN in maize. Some tassel architecture-related genes were cloned based on maize mutants previously. For example, RAMOSA2 (RA2), and RAMOSA3 (RA3), which function in the ramosa pathway, regulate the identity and determinacy of inflorescence axillary meristems. Additionally, those ramosa mutants show an increase in tassel branch number [13,14]. The two fasciated ear mutants, fasciated ear2 (fea2) and fasciated ear3 (fea3), also showed a larger tassel branch number. FEA2 and FEA3 encode membrane-localized leucine-rich repeat receptors and involve in the CLAVATA pathway [15,16]. UNBRANCHED2 (UB2), UB3, and TSH4 all encode SBP-box genes and function as redundant factors that limit the rate of cell differentiation to the lateral primordia. Loss-of-function mutations in ub2, ub3, and tsh4 result in smaller TBN [17,18]. Barren inflorescence2 (bif2) mutation has a reduced number of branches and spikelets in tassels. BIF2 encodes a serine-threonine protein kinase, which is modified the transition from inflorescence meristem (IM) to spikelet pair meristem (SPM) and indirectly affected the development of branch meristems [19,20]. However, because of negative pleiotropy, these genes cloned by mutants are difficult to use in maize breeding. Recently, a natural variation related to tassel branch number was identified in the 5 -UTR (untranslated region) of Q Dtbn1 . Q Dtbn1 encodes a Kelch repeat-containing F-box protein, negatively regulates TBN with a dominant model, and participates in the ABA signaling pathway [21]. Despite these studies, there are still very few genes cloned through natural variation.
In this study, a B73 × CML247 RILs population was used to detect QTLs for TBN. Four QTLs were detected and colocated with the TBN-related associated loci that were reported by Wu et al. (2016), using two NAM populations to identify candidate genes [12]. ZmPAT7 was identified as a major candidate gene through bioinformatics and expression analysis. ZmPAT7 negatively regulated TBN, and its knockout lines had more TBN. A key SNP variation in the promoter of ZmPAT7 may confer the expression of ZmPAT7 and shows significant association with TBN. Moreover, 17 distant eQTLs associated with the expression of ZmPAT7 were identified in expression quantitative trait loci (eQTL) analysis.

QTL Mapping and Identification of Candidate Genes for TBN
The two inbred lines, B73 and CML247, display great differences in tassel branch number. CML247 plants had more TBN than B73 plants, and the B73 × CML247 RIL population also revealed considerable variations in TBN ( Figure S1a). Therefore, we used the RIL population to conduct subsequent QTL mapping. A total of four QTLs were detected above a cutoff of LOD > 2.5, including two QTL (qTBN2a and qTBN2b), one QTL (qTBN4), and one QTL (qTBN6) on chromosomes 2, 4, and 6, respectively ( Figure 1, Table 1). Individual QTLs explained phenotypic variation ranging from 5.21% to 16.22%, and qTBN4 showed the largest effect among the four QTLs (Table 1). The number of genes in the interval of QTLs varied from 217 to 2510 (Table 1). We further checked whether the known cloned genes regulating the development of male inflorescence were located in QTL regions. We found that FASCIATED EAR2 (FEA2), involved in the CLAVATA pathway to control the size of the inflorescence meristem, was colocated with qTBN4. zfl2, whose mutant exhibited disruptions of floral organ identity and patterning, and defects of inflorescence architecture, and in the vegetative to reproductive phase transition [22], they were colocated with qTBN2a.
control the size of the inflorescence meristem, was colocated with qTBN4. zfl2, whose mutant exhibited disruptions of floral organ identity and patterning, and defects of inflorescence architecture, and in the vegetative to reproductive phase transition [22], they were colocated with qTBN2a.   [12]. It was found that 25 among the 549 QTNs were colocated in the four QTLs detected in this study, including 15 QTNs colocated with qTBN2a, 4 QTNs colocated with qTBN2b, 5 QTNs colocated with qTBN4, and 1 QTN colocated with qTBN6. Subsequently, we identified a total of 36 candidate genes underlying TBN (Table 2). Furthermore, we compared the expression of those candidate genes in shoot apical meristem (SAM) between the two parents, i.e., B73 and CML247 ( Figure S1b). Six candidate genes (GRMZM2G302701, GRMZM5G897776, GRMZM2G162266, GRMZM2G163067, GRMZM2G010011, and GRMZM2G423636) showed high expression and also differential expression between B73 and CML247, which might play a role in regulating maize TBN. The three candidate genes, i.e., GRMZM2G162266, GRMZM2G163067, and GRMZM2G423636, encode uncharacterized proteins. GRMZM2G010011 encodes an Sacyltransferase. GRMZM2G302701 encodes a nucleolar preribosomal-associated protein 1. GRMZM5G897776 encodes a starch synthase 3 involved in starch biosynthesis.   Previously, Wu et al. (2016) performed a genome-wide association study using the CN-NAM and US-NAM populations for tassel architecture and identified a total of 549 QTNs for TBN [12]. It was found that 25 among the 549 QTNs were colocated in the four QTLs detected in this study, including 15 QTNs colocated with qTBN2a, 4 QTNs colocated with qTBN2b, 5 QTNs colocated with qTBN4, and 1 QTN colocated with qTBN6. Subsequently, we identified a total of 36 candidate genes underlying TBN (Table 2). Furthermore, we compared the expression of those candidate genes in shoot apical meristem (SAM) between the two parents, i.e., B73 and CML247 ( Figure S1b). Six candidate genes (GRMZM2G302701, GRMZM5G897776, GRMZM2G162266, GRMZM2G163067, GRMZM2G010011, and GR-MZM2G423636) showed high expression and also differential expression between B73 and CML247, which might play a role in regulating maize TBN. The three candidate genes, i.e., GRMZM2G162266, GRMZM2G163067, and GRMZM2G423636, encode uncharacterized proteins. GRMZM2G010011 encodes an S-acyltransferase. GRMZM2G302701 encodes a nucleolar preribosomal-associated protein 1. GRMZM5G897776 encodes a starch synthase 3 involved in starch biosynthesis.

ZmPAT7 Negatively Controls Tassel Branch Number in Maize
Previous studies suggested that S-acyltransferases play an important role in regulating plant architecture in plants [23][24][25][26]. Therefore, we focused on GRMZM2G010011 to further validate its function. Previously, 38 Protein S-acyltransferases (PAT) genes have been identified in the maize genome, and GRMZM2G010011 was named ZmPAT7 [27]. ZmPAT7 harbors a conserved Asp-His-His-Cys (DHHC) domain that is representative of S-acyltransferases and has four transmembrane motifs ( Figure S2a,b). Phylogenetic analyses of ZmPAT7 with high homology protein from other plants showed that ZmPAT7 was closely related with PAT proteins in monocot ( Figure 2a). To confirm the function of ZmPAT7 underlying TBN, we obtained two knockout lines (KO#1 and KO#2) using CRISPR/Cas9 technology. Gene editing of ZmPAT7 was screened by PCR amplification and Sanger sequencing of the target region. KO#1 carrying a 1 bp deletion and KO#2 carrying a 27 bp deletion were used for subsequent analysis (Figure 2b). Compared with the wild-type line, the two loss-of-function lines of ZmPAT7 showed a significant increase in TBN (Figure 2c,d).

Natural Alleles of ZmPAT7
To identify key variations in ZmPAT7 associated with TBN, we obtained the polymorphic variants in a~12.5 kb genomic region, covering the promoter and gene body of ZmPAT7, from 1604 resequenced maize inbred lines. We identified 200 polymorphic variants (SNPs and indels) ( Figure 3a) and found that three SNPs (−101, 3771, and 8855 bp) were significantly associated with TBN. The three SNPs were located in the promoter, the fifth intron region, and downstream of the gene body, respectively. Based on the three significant variants, we classified the 1604 maize genotypes into 4 haplotype groups (Figure 3b). The lines with Hap1 or Hap2 have lower TBN than those with Hap3 or Hap4, and no significant differences were observed between lines with Hap1 and Hap2 or between lines with Hap3 and Hap4. Therefore, we designated Hap1 or Hap2 as a favorable haplotype, which is consistent with the haplotype groups classified only by SNP (−101). The SNP (−101) may be a key variant in regulating TBN. For example, the inbred lines Zheng58 and Chang7-2, two parents of hybrid Zhengdan958 widely planted in China, had significant differences in TBN (Figure 3c,d). Zheng58 with Hap1 had a smaller TBN (5.5 ± 1.1) than Chang7-2 with Hap4 (18 ± 1.9). Furthermore, we chose 20 lines with Hap1 and 20 lines with Hap4 to conduct gene expression analysis. ZmPAT7 exhibited higher expression in low TBN lines with Hap1 than that in high TBN lines with Hap4 ( Figure 3e). These results suggested that the expression of ZmPAT7 may contribute to the phenotypic variation in maize tassel branch number.

Natural Alleles of ZmPAT7
To identify key variations in ZmPAT7 associated with TBN, we obtained the polymorphic variants in a ~12.5 kb genomic region, covering the promoter and gene body of ZmPAT7, from 1604 resequenced maize inbred lines. We identified 200 polymorphic variants (SNPs and indels) ( Figure 3a) and found that three SNPs (−101, 3771, and 8855 bp) were significantly associated with TBN. The three SNPs were located in the promoter, the fifth intron region, and downstream of the gene body, respectively. Based on the three significant variants, we classified the 1604 maize genotypes into 4 haplotype groups (Figure 3b). The lines with Hap1 or Hap2 have lower TBN than those with Hap3 or Hap4, and no significant differences were observed between lines with Hap1 and Hap2 or between lines with Hap3 and Hap4. Therefore, we designated Hap1 or Hap2 as a favorable haplotype, which is consistent with the haplotype groups classified only by SNP (−101). The SNP (−101) may be a key variant in regulating TBN. For example, the inbred lines Zheng58 and Chang7-2, two parents of hybrid Zhengdan958 widely planted in China, had significant differences in TBN (Figure 3c,d). Zheng58 with Hap1 had a smaller TBN (5.5 ± 1.1) than Chang7-2 with Hap4 (18 ± 1.9). Furthermore, we chose 20 lines with Hap1 and 20 lines with Hap4 to conduct gene expression analysis. ZmPAT7 exhibited higher expression in low TBN lines with Hap1 than that in high TBN lines with Hap4 (Figure 3e). These results suggested that the expression of ZmPAT7 may contribute to the phenotypic variation in maize tassel branch number.

Expression Quantitative Trait Loci Analysis of ZmPAT7
We used 223 inbred lines with large phenotypic variations in TBN (from 1.97 to 23.37) (Table S1, Figure S3) to analyze the expression level of ZmPAT7. The expression of ZmPAT7 showed a considerable variation and negatively correlated with TBN (r = −0.29, p = 2.7 × 10 −4 ) (Figure 4a). Association analysis was conducted using the expression of ZmPAT7 as phenotypic value and a total of 49 significantly associated SNPs were identified (Figure 4b). Based on the significant SNPs, a total of 17 eQTLs associated with the expression of ZmPAT7 were identified. These included nine on chromosome 2, three on chromosome 3, three on chromosome 5, one on chromosome 6, and one on chromosome 7 (Table 3). All the 17 eQTLs were defined as "distant eQTLs" following the previous method proposed by Fu et al. (2013) and   [28,29]. In total, 44 genes were found in LD distance (200 kb) adjacent to the SNPs (Table 3), among which some transcription factor genes, e.g., Zm00001d038207 (ZmNAC3), may be involved in regulating ZmPAT7. Interestingly, we also found an SNP (−274, T/A) in the promoter of ZmPAT7 between B73 and CML247, which results in the cis-element NAC motif change. The NAC motif (TCTTGACC) existed in B73, while the NAC core promoter element was absent in CML247 ( Figure S4). These results suggested that ZmNAC3 may be a major factor that mediated the expression of ZmPAT7.

Expression Quantitative Trait Loci Analysis of ZmPAT7
We used 223 inbred lines with large phenotypic variations in TBN (from 1.97 to 23.37) (Table S1, Figure S3) to analyze the expression level of ZmPAT7. The expression of ZmPAT7 showed a considerable variation and negatively correlated with TBN (r = −0.29, p = 2.7 × 10 −4 ) (Figure 4a). Association analysis was conducted using the expression of ZmPAT7 as phenotypic value and a total of 49 significantly associated SNPs were identified (Figure 4b). Based on the significant SNPs, a total of 17 eQTLs associated with the expression of ZmPAT7 were identified. These included nine on chromosome 2, three on chromosome 3, three on chromosome 5, one on chromosome 6, and one on chromosome 7 (Table 3). All the 17 eQTLs were defined as "distant eQTLs" following the previous method proposed by Fu [28,29]. In total, 44 genes were found in LD distance (200 kb) adjacent to the SNPs (Table 3), among which some transcription factor genes, e.g., Zm00001d038207 (ZmNAC3), may be involved in regulating ZmPAT7. Interestingly, we also found an SNP (−274, T/A) in the promoter of ZmPAT7 between B73 and CML247, which results in the cis-element NAC motif change. The NAC motif (TCTTGACC) existed in B73, while the NAC core promoter element was absent in CML247 ( Figure S4). These results suggested that ZmNAC3 may be a major factor that mediated the expression of ZmPAT7.

Discussion
In this study, we mapped four QTLs for TBN by using a high-quality recombination map of the B73 × CML247 RIL population. In previous studies, all identified QTLs were colocated with the QTLs/QTNs detected in other populations [6,9,11,12]. qTBN2a mapped at 7.85-14.72 Mb interval on chromosome 2 in our study was also reported previously. For example, Q32 CN-NAM and Q46 US-NAM identified by Wu et al. (2016) using joint linkage analysis for TBN [12] were colocated with qTBN2a. Zhao et al. (2017) also identified an mQTL2-1 on chromosome 2 using meta-QTL analysis [6], which overlapped with qTBN2a in our study. In addition, a tassel dry weight (TW) related QTL, qTW2_1, was also colocated with qTBN2a [9]. Q50 US-NAM and four QTNs reported by Wu et al. (2016) [12] were colocated with qTBN2b, and Q44 CN-NAM reported by the same study was colocated with qTBN6. Eight QTNs and five QTNs reported by Brown et al. (2011) [11] and Wu et al. (2016) [12], respectively, fell in the interval of qTBN4. Furthermore, some maize inflorescence developmental genes were located in those QTLs detected in this study. For example, zfl2, regulating the development of inflorescence architecture, was colocated with qTBN2a. In addition, fea2, controlling the size of the inflorescence meristem, was colocated with qTBN4. These comparisons suggested that the QTLs detected in the present study were important genomic regions controlling TBN.
A large number of QTLs for TBN have been identified by using bi-parental populations [6,8,9,[30][31][32][33]. These QTLs have generally large genomic intervals, which is difficult to find candidate genes quickly. The interval of QTLs in our study is also large, containing many genes (Table 1). On another front, many QTNs for TBN were identified using GWAS [11,12,34]. Combining GWAS and linkage mapping would be beneficial to quickly identify candidate genes.  [12]. A total of 25 QTNs were colocated with the 4 QTLs obtained in this study, and 36 candidate genes were then identified.
Transcriptome data can be used in narrowing the range of candidate genes for specific traits. Thus, we used the expression data in B73 and CML247 shoot apical meristem (SAM) from https://qteller.maizegdb.org/ (accessed on 18 May 2020) [35] for identifying the most possible candidate genes. Six candidate genes (GRMZM2G302701, GRMZM5G897776, GRMZM2G162266, GRMZM2G163067, GRMZM2G010011, and GR-MZM2G423636) were identified to further explore the regulation mechanisms of TBN, among which GRMZM2G010011 (ZmPAT7) encodes an S-acyltransferase.
S-acylation is a reversible post-translational modification process occurring in eukaryotes that regulates tracking, regulation, signaling, membrane association, and target protein functions [36]. Several studies have reported that S-acyltransferases encoded by the PAT gene family play an important role in regulating plant architecture and other traits in plants. For instance, OsPAT15 regulates plant architecture by altering the tiller in rice [24], and overexpression of OsPAT15 results in increased branch and seed yield in Brassica napus [37]. Knockout lines of PbPAT14 exhibited the dwarf yellowing phenotype, and the function of PbPAT14 was related to the ABA pathway [25]. In addition, Hemsley et al. (2005) reported an S-acyltransferase gene in Arabidopsis, Tip Growth Defective 1 (TIP1)/PAT24, which exhibits impairment in root hair elongation and defective growth of pollen tubes [23]. Zhang et al. (2020) reported that ZmTIP1, which encodes a functional S-acyltransferase, plays a positive role in regulating the length of root hairs and drought tolerance in maize [26]. Two homologous proteins, PAT13 and PAT14, are cooperatively involved in regulating precocious leaf senescence in Arabidopsis [38,39]. PAT10 is localized on Golgi and tonoplast membranes, and those mutations performed pleiotropic defects in plant development [40,41]. PAT4 was also reported to mediate root hair growth in Arabidopsis [42]. AtPAT21, another PAT reported recently, participates in both male and female gametogenesis [43]. Although ZmPAT7 harbors a conserved DHHC domain that is representative of S-acyltransferases (Supplementary Figure S2a,b), and phylogenetic analyses of ZmPAT7 showed that it was highly conserved in monocots, functions of these S-acyltransferase proteins that are close to ZmPAT7 are still unknown.
In our study, the two loss-of-function lines of ZmPAT7 showed significant increases in TBN (Figure 2c,d). ZmPAT7 exhibited higher expression in low TBN lines than that in high TBN lines (Figure 3d), suggesting that the expression of ZmPAT7 negatively regulated TBN in maize. We identified a key variant SNP (−101) in the ZmPAT7 promoter which caused the change in expression level and was further verified in the inbred lines of Zheng58 and Chang7-2 (Figure 3b,c). Additionally, the sequences of ZmPAT7 promoter in B73 and CML247 also showed several variants (Supplementary Figure S4), e.g., the SNP (−157, A/C) resulting in the ERF element difference, the SNP(−274, T/A) in the cis-element NAC motif change, and the SNP(−433, G/T) in the MYB element difference. Taken together, we speculated that ZmPAT7 plays an important role in natural variations in TBN in maize.

Plant Materials
The B73 × CML247 RILs population containing 193 lines was used for QTL mapping. A large association population containing 1604 maize inbred lines comprising 1227 from China, 344 from the United States, and 33 from other countries [44] was used for candidate gene association analysis. Furthermore, we chose 223 inbred lines that showed considerable variations in TBN from the association panel to conduct expression assay and eQTL analysis. To detect the expression of ZmPAT7, we grew the 223 inbred lines in a field in Beijing, China, in the summer of 2019, and sampled the leaf at the V5 stage.

QTL Mapping and Identification of Candidate Genes
The phenotypic data of TBN for the US-NAM RILs were collected in eight environments as described by Brown et al. (2011) [11], which were transformed using the Box-Cox function in R software, and used to calculate best-unbiased predictor (BLUP). The BLUP values of TBN in B73 × CML247 RILs were used for QTL mapping. A high-quality recombination map containing 2183 bin markers of the B73 × CML247 RILs was reported by Li et al. (2015) [45]. The total linkage distance of the genetic map is 1708.44 cM, and the average interval between adjacent bins was 0.78 cM [45]. A total of 193 RILs were matched between the phenotype and genotype. QTL analysis was conducted using the WinQTLCart cartographer v2.5 software (North Carolina State University, Raleigh, NC, USA) with composite interval mapping (CIM) method [46], with a 1000 permutation test at 95% confidence level to determine the optimal log of odds (LOD) threshold values. Then, the 2-LOD drop method was applied for defining the QTL confidence interval. Then, we identified the colocated loci with the QTNs reported by Wu et al. (2016) [12]. Genes falling within 50 kb of colocated loci were considered candidate genes. Furthermore, the expression data in B73 and CML247 shoot apical meristem (SAM) of candidate genes was extracted from https://qteller.maizegdb.org/ (accessed on 18 May 2020) [35] and used to find differential expression genes between B73 and CML247. The genes with high expression in SAM and the ratio of expression between B73 and CML247 greater than 1.5 (B73/CML247 > 1.5 or CML247/B73 > 1.5) were considered candidate genes.

Identification of the CRISPR/Cas9 Mutants of ZmPAT7
The CRISPR/Cas9 mutants of ZmPAT7 were ordered from http://www.wimibio.com/ maize.html (accessed on 5 May 2019) and the methods of vector construction, genetic transformation for creating CRISPR/Cas9 mutants were reported by Liu et al. (2020) [47]. Briefly, the single guide RNA (sgRNA) sequence is AGGGCTCAGAGTGAACCTCC. The vector pCPB-ZmUbi-hspCas9 was first linearized by HindIII. Separately, ZmU6 and the sgRNA scaffold of ZmPAT7 were amplified through overlapping PCR with a homologous arm. Additionally, homologous arms that match linearized pCPB-ZmUbi-hspCas9 were also added to the insertion fragment in the overlap PCR. Finally, the gene-specific insertion fragments were incorporated into the pCPB-ZmUbi-hspCas9 vector. Additionally, the maize inbred line KN5585 was then used for Agrobacterium tumefaciens-mediated transformation of immature embryos. We planted the T1 lines in the field to identify the homozygous positive lines. Genomic editing of ZmPAT7 was screened by PCR amplification and Sanger sequencing of the target region. The PCR primers were cr-F: CAATCTTGACAAAGAAACC and cr-R: TTGAAATGATAGAACGCAC. The T2 homozygous positive lines and its wild-type line KN5585 were phenotyped at Sanya (18 • N, 109 • E) in the winter of 2020. The experiment in the field followed a randomized block design, with three replicates.

Candidate Gene-Based Association Analysis
The phenotypic data of TBN for the 1604 inbred lines were identified at five locations in 2018 and 2019, which were used to calculate the best-unbiased predictor (BLUP). The polymorphic variants in a~12.5 kb region covering 2.5 kb upstream, 9.5 kb gene body of ZmPAT7 and 0.5 kb downstream region were extracted from 1604 inbred lines resequenced data [44]. Those variants with minor allele frequency (MAF) over 5% were identified and used for candidate gene-based association analysis. Association analysis was performed using a mixed linear model (MLM) in Tassel v3.0 (Cornell University, Ithaca, NY, USA), with the population structure and kinship analyzed in Tassel v5.0 [48]. A Bonferroni-corrected significant association threshold (p ≤ 0.1/200 = 5.0 × 10 −4 ) was used to determine the significantly associated variants.

Expression Quantitative Trait Loci Analysis
The leaves at V5 of 223 inbred lines were isolated for qRT-PCR. A total of 8,974,340 SNPs was identified with a minor allele frequency (MAF) > 0.05. The Kinship (K) matrix of pairwise genetic distances was calculated by EMMAX and was used as the variancecovariance matrix of the random effects for association analysis. The association analysis using a mixed linear model (MLM) was performed with the EMMAX software package. The suggestive threshold was 1.0 × 10 −5 , which was used to identify significantly associated signals. The significant SNPs were defined as candidate eQTLs, and if there were more than two SNPs in a 5 kb region, the SNPs were regarded as a cluster, and the most significant SNPs in the clusters were considered as candidate eQTLs. If two candidate eQTLs of a single gene were under LD (r 2 > 0.1), the less significant eQTL was removed. The eQTLs were divided into two types using 20 kb as a cutoff to define "local eQTLs" versus "distant eQTLs". The methods mentioned above were similar to those in previous studies [28,29].

RNA Isolation and Quantitative Reverse Transcription PCR (qRT-PCR)
Total RNAs were extracted from plant tissues using the Plant Total RNA Isolation Kit (Genebetter, Beijing, China). First-strand cDNA was synthesized using the Uni One-Step gDNA Removal and cDNA Synthesis SuperMix Kit (TranScript, Beijing, China). Quantitative real-time PCR (qRT-PCR) was performed using the Taq Pro Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) with a QuanStudio 3 Real-Time PCR System cycler (Applied Biosystems, San Francisco, CA, USA). The maize GAPDH gene was used as the internal control, and the relative expression was calculated using the 2 −∆∆Ct method. The PCR conditions consisted of an initial denaturation step at 95 • C for 30 s, followed by 40 cycles at 95 • C for 10 s, and 60 • C for 30 s. The primers of ZmPAT7 and GAPDH were qPAT7 (qPAT7-f: CCAGTATGAGGCAGACAGTGA, qPAT7-R: TGGATGGTGCTCGGCTATG) and qGAPDH (qGAPDH-F: CCCTTCATCACCACGGACTAC qGAPDH-R: AACCTTCTTG-GCACCACCCT), respectively.