A common wheat (Triticum aestivum) variety, the Chinese Spring, was used for cloning of TaPYL4 gene sequences. Thirty-two wheat accessions (Supplementary Table 1) with wide variation screened by 209 SSR markers (Zhang et al., 2013) were used to get gene sequences for detecting SNPs of TaPYL4s.

Three wheat populations were employed in this research. Population 1 (323 accessions) (Wang et al., 2019a) was used for agronomic character analysis. This population consisting of 275 modern cultivars, 36 advanced lines, and 12 landraces was also employed for association analysis. Population 2 and Population 3 consisted of 157 landraces and 348 modern cultivars (Zhang et al., 2002), respectively. The latter two populations were mainly from a Chinese wheat microcore germplasm, accounting for more than 70% of the genetic diversity of the whole Chinese wheat germplasm resources (Hao et al., 2011a). Haplotype frequencies among 10 wheat-producing areas in China were determined using Population 2 and Population 3 (Zhang et al., 2002).

Population 1 was planted in 2015–2017 in Shunyi (40°23′N; 116°56′E) and Changping (40°13′N; 116°13′E) field experiment stations of the Institute of Crop Sciences, Chinese Academy of Agricultural Sciences, Beijing. Field water management was divided into two conditions: rain-fed (drought stressed, DS) and well-watered (WW). DS referred to that no artificial irrigation was given in the whole growing season of wheat, and the growth was completely fed by rain. WW was designed to irrigate the field with 750 m³ha^–1 (75 mm) water at pre-wintering, jointing, flowering, and grain-filling stages, respectively. Moreover, a heat stress experiment (HS) was applied to wheat plants of 1-week post-anthesis by placing a plastic film supported by steel frames over the plots in the field at Shunyi. Wheat phenotypic data under 16 environments (E1-E16) were obtained. E1-E16 represented the environments at Shunyi in 2015 under DS + HS, DS, WW + HS and WW, at Shunyi in 2016 under DS + HS, DS, WW + HS and WW, at Changping in 2016 under WW and DS, at Shunyi in 2017 under DS + HS, DS, WW + HS and WW, at Changping in 2017 under DS and WW (Li et al., 2019a), respectively (Supplementary Table 2).

The agronomic traits including plant height (PH), spike length (SL), peduncle length (PLE), length of penultimate node (LPN), number of spikes per plant (NSP), number of spikelet per spike (NSS), number of grains per spike (NGS), and 1,000-grain weight (TGW) were measured under 16 environments. Population 2 and 3 were grown at Luoyang (36°61′N; 112°45′E) in Henan Province in 2002 and 2005, and also grown at Shunyi (40°23′N; 116°56′E) in Beijing in 2010 (Zhang et al., 2017). In all cases, each variety was planted on a 2-m, 4-row plot with 30 cm row spacing of 40 seeds in each row. Phenotypic sampling was performed for 5 plants grown in the middle of each plot at maturity.

The reference sequence (Accession number MG273654 in the NCBI database) of TaPYL4 was obtained according to Mega’s research (Mega et al., 2019). Genome-specific primers 2A-F1/R1, 2B-F1/R1, and 2D-F1/R1 were designed to amplify the promoter and coding region sequences of Chinese Spring wheat using Primer Premier 5.0 software (Supplementary Table 3). MEGA7 software was used to construct a neighbor joining phylogenetic tree.

All genomic DNA used in this study was extracted from the wheat by CTAB method (Stewart and Via, 1993). Gene promoter and coding region sequences of TaPYL4s were amplified from the 32 diverse wheat accessions (Li et al., 2019b), and subsequently ligated into pEASY-Blunt vector, respectively (Trans Gen Biotech, Beijing). Positive clones were randomly selected, and the plasmid from each clone was extracted for sequencing.

According to the sequence information of TaPYL4-2A, TaPYL4-2B, and TaPYL4-2D, genome-specific sequencing primers CX-A-R1/R2, CX-B-R1/R2, and CX-D-F1/R1 (Supplementary Table 3) were designed by Primer 5 and using a 3730 XI DNA Analyzer (ABI) sequencing. Seqman software in DNASTAR Lasergene 7.1 package was used for sequence splice and assembly. Meg Align software was used for sequence polymorphism analysis.

According to the two polymorphic sites detected, two molecular markers (TaPYL4-2A-dCAPS and TaPYL4-2B-dCAPS) were established at A genome −1635 bp (G/A) and B genome −1146 bp (G/C) upstream of the ATG translation start site, respectively. These dCAPS primers were designed with an available program dCAPS Finder 2.0. The primers were named 2A-SalI-F/R and 2B-BamHI-F/R, respectively (Supplementary Table 3). The first PCR product obtained with A genome-specific primers and B genome-specific primers, respectively, were served as the template for the second round of PCR. The annealing temperature and elongation time were determined by primer pairs and lengths of the expected PCR products. The PCR products were digested with the appropriate restriction enzyme at 37°C for 1.5 h and further separated on 4% agarose gel electrophoresis.

Genotype scanning was performed on the Population 1 using the molecular markers (TaPYL4-2A-dCAPS and TaPYL4-2B-dCAPS) developed. Population 1 was performed on Wheat 660 K SNP Array, which consisted of 630,517 SNPs (Li et al., 2019b). By removing nucleotide variations with missing rates ≥ 0.2 and minor allele frequency (<0.05), 395,681 SNPs were eventually used to detect the structure of population 1 by software STRUCTURE 2.3.4 (Li et al., 2019a). Using the general liner model (GLM) of TASSEL (Version 5.2.51) software and population structure (Q) as covariable, correlation analysis was conducted for phenotypic traits and genotypes of TaPYL4-2A and TaPYL4-2B. Associations were considered as significant at P < 0.05. Analysis of variance by SPSS (version 19.0) software was used to detect different effects of haplotype on the agronomic traits.

Based on the above sequencing results, analysis of cis-acting element in Plant CARE¹ was employed to predict the cis-acting element in the promoter region of around 2,000 bp upstream of TaPYL4-2A and TaPYL4-2B genes, respectively.

The expression pattern of TaPYL4-2A was analyzed using the wheat accession Chinese Spring. The spatiotemporal expression patterns of TaPYL4-2A were examined for different tissues (shoot and root at seedling stage, flag leaves, spikes, penultimate nodes, root bases, and root at heading stage). Twelve accessions of two different haplotypes were randomly selected from Population 1 (Supplementary Table 4). The gene expression levels were evaluated in shoot (St) for 20 d-old seedlings. Three biological replications were performed with each repeat at least three technical duplications.

Yeast one-hybrid assays were employed to verify the binding of TaARF4 to TaPYL4-2A in different haplotype promoter regions. NotI was used to digest the pB42AD vector to obtain TaARF4 fragment with a restriction enzyme site. Next, the TaARF4 was cloned into the vector with seamless clonal enzyme. Two fragments of different haplotype promoter regions of TaPYL4-2A (−1890 to −1440 bp) were cloned into pLacZi as the reporter gene plasmid, respectively. The plasmid of pB42AD carrier with TaARF4 was co-transformed into EGY48 in yeast with 2 reporters by standard yeast transformation method, respectively (Lin et al., 2007). The yeast was grown on SD/-Trp/-Ura medium with x-gal. LacZ staining was used to detect the interaction between TaARF4 and promoter fragments (Wang et al., 2016). This experiment was conducted in triplicate.

Primers TaARF4-1,300-F/R (Supplementary Table 3) were used to amplify TaARF4 full-length cDNA. The cDNA was further cloned into effector vector pCAMBIA1300 driven by CaMV35S. TaPYL4-2A promoter fragments of different haplotypes were amplified by primers PYL4-LUC-F1/R1 (Supplementary Table 3) and then cloned into the reporting vector pGreenII0800-LUC. The plasmid of pCAMBIA1300 vector with TaARF4 and different reporter genes were separately transformed into GV3101 cells containing pSoup. The Agrobactium GV3101 containing the recombinant vector was used to transform tobacco leaves of 4 weeks by co-osmosis. The luciferase activity of firefly luciferase (LUC) and Renilla luciferase (REN) were determined at 48 h after osmosis using a Dual Glo^® Luciferase Assay System (E2920, Promega) with a multimode reader (TriStar² S LB942). The ratio of LUC to REN was used to calculate promoter activity.

Three sub-genomic sequences of TaPYL4 gene were isolated from a wheat accession, the Chinese Spring, namely as TaPYL4-2A, TaPYL4-2B, and TaPYL4-2D, respectively. The gene coding regions from A, B, and D sub-genomes are 654, 654, and 645 bp in length, respectively. Like other PYL family members, TaPYL4 proteins were highly conserved (Supplementary Figure 1A). The phylogenetic tree showed that the TaPYL4s were more similar to OsPYL4 (Supplementary Figure 1B).

In the present study, SNPs in the coding and promoter regions of TaPYL4-2A, -2B, and -2D genes were detected by sequencing 32 diverse wheat accessions. No SNP was identified in the coding regions of TaPYL4-2A and TaPYL4-2B. However, nine SNPs were observed in TaPYL4-2A promoter region while two SNPs existed in TaPYL4-2B promoter region. For TaPYL4-2D, neither coding sequence nor promoter region was detected to have SNP. Therefore, no further investigation was on TaPYL4-2D in this study. According to the SNPs in TaPYL4-2A (Figure 1A) and TaPYL4-2B (Figure 2A) promoter regions, two haplotypes were identified for TaPYL4-2A (Hap-2A-1 and Hap-2A-2) and TaPYL4-2B (Hap-2B-1 and Hap-2B-2), respectively.

To accurately distinguish two haplotypes of TaPYL4-2A, a molecular marker (dCAPS) was developed based on the SNP (G/A) at −1635 bp of the gene, namely as TaPYL4-2A-dCAPS. The marker included two mismatches (AC/TT) in the downstream primer that generated the recognition site for the restriction enzyme SalI (Figure 1B). In a similar way, a dCAPS marker was also developed with restriction enzyme BamHI based on TaPYL4-2B, namely as TaPYL4-2B-dCAPS (Figure 2B).

Association analysis was conducted on agronomic characters of Population 1 under 16 environments in three growth seasons. For TaPYL4-2A, Hap-2A-1 accounted for 48.9% in the Population 1 (Supplementary Figure 2). The association analysis of TaPYL4-2A haplotypes and agronomic traits revealed that TaPYL4-2A was significantly associated with SL, PLE, and PH (Table 1). The statistical analysis showed that the SL of Hap-2A-1 was significantly shorter than that of Hap-2A-2 under 16 environmental conditions (Figure 3A). Under 14 environmental conditions, the PLE of Hap-2A-1 was significantly shorter than that of Hap-2A-2 (Figure 3B). Under eight environmental conditions, the PH of Hap-2A-1 was significantly shorter than that of Hap-2A-2 (Figure 3C). In summary, the values measured for plant growth related traits of Hap-2A-1 were lower than that of Hap-2A-2.

For TaPYL4-2B, Hap-2B-1, and Hap-2B-2 accounted for 75.5 and 24.5% in Population 1 (Supplementary Figure 2), respectively. Significant association between TaPYL4-2B haplotypes and PH was identified under 13 environments (Table 2). The PH of Hap-2B-1 was significantly shorter than that of Hap-2B-2 (Figure 4A). Similarly, TaPYL4-2B haplotypes were associated with PLE. Hap-2B-1 had significantly shorter PLE compared with Hap-2B-2 in 14 environments (Figure 4B). According to the description above, the values measured for plant growth related traits of Hap-2B-1 were lower than that of Hap-2B-2.

Since both TaPYL4-2A and TaPYL4-2B were associated with plant growth related traits (PH and PLE), we speculated that the TaPYL4-2A and TaPYL4-2B haplotypes may have a combinational effect on plant growth related traits. The present association studies showed that both TaPYL4-2A and TaPYL4-2B haplotypes can be divided into low PH and PLE haplotypes (Hap-2A-1 and Hap-2B-1) and high PH and PLE haplotypes (Hap-2A-2 and Hap-2B-2). Wheat accessions presenting these four combinations were divided into groups of low-/ low-, low-/ high-, high-/ low-, and high-/high-PH and PLE haplotypes, which were named as Hap-AB1, Hap-AB2, Hap-AB3, and Hap-AB4 (Supplementary Table 5), respectively. Among the four combinations, Hap-AB1 was associated with the lowest PLE (Figure 5A) and PH (Figure 5B) in Population 1, while Hap-AB4 associated with the highest PH and PLE. The plant growth related traits of Hap-AB2 was lower than that of Hap-AB3. These results demonstrated the TaPYL4-2A haplotypes had a larger contribution to plant growth related traits than the TaPYL-2B haplotypes’ effects. Furthermore, an additive effect of the TaPYL4-2A and TaPYL4-2B haplotypes was detected to work on the regulation of plant growth.

Cis-acting elements in promoter region could reflect the function and regulatory mechanism of the genes. In order to further clarify the influence of polymorphic sites located in TaPYL4-2A and TaPYL4-2B promoter region, we used Plant CARE² to analyze the cis-acting elements in promoters of two haplotypes detected for TaPYL4-2A and TaPYL4-2B, respectively. Interestingly, an ARF binding site, TGA-element (AACGAC), was detected at −1635 bp upstream of the ATG translation starting point of Hap-2A-1. This cis-acting element was absent in the promoter of Hap-2A-2 due to a nucleotide variation (Supplementary Figure 3). However, in the promoter regions of two haplotypes of TaPYL4-2B, there was no difference in cis-acting elements related to plant growth and development (Supplementary Figure 4).

Since the nucleotide variation site is located in the TaPYL4-2A promoter region, we further examined the impacts of this variation on the gene expression. The expression pattern of TaPYL4-2A gene was identified by real-time PCR at different stages of growth, showing that TaPYL4-2A was constitutively expressed in various tissues tested (Figure 6A).

The 12 wheat accessions representing two haplotypes of TaPYL4-2A were randomly selected from Population 1 and planted in the field. After 20 days of growth, shoot samples were taken for measuring the expression of TaPYL4-2A. Notably, the expression level of TaPYL4-2A in Hap-2A-1 wheat accessions was higher than that in Hap-2A-2 wheat accessions (Figure 6B).

Association analysis and expression patterns indicated that the expression level of TaPYL4-2A was high in Hap-2A-1 wheat accessions with shorter SL and PLE as well as shorter PH. On the contrary, the low expression of Hap-2A-2 was detected in the accessions with longer SL and PLE as well as taller PH. Therefore, TaPYL4-2A gene expression might be negatively correlated with plant growth-related traits.

Previous reports showed that ARF7 and ARF19 could be bound to the TGA-element in Arabidopsis (Huang et al., 2018). TaARF4 was related to plant height in wheat (Wang et al., 2019b). Therefore, we speculated that TaARF4 may regulate the expression of TaPYL4-2A. Yeast one-hybrid assays were carried out to determine whether TaARF4 binds to the promoter region of TaPYL4-2A. The TaARF4 DNA binding region was cloned into pB42AD vector. The promoter regions of the haplotype TaPYL4-2A with and without the TGA-element were cloned into the pLacZi, separately (Figure 7A). In order to identify the interaction between promoters and TaARF4, the plasmid containing each of these 2 promoter regions was transformed into yeast EGY48 with TaARF4 plasmid, respectively. Positive yeast EGY48 cells were cultured on the selective medium of x-gal without Ura and Trp. In the presence of TGA-element in Hap-2A-1 of TaPYL4-2A, TaARF4 were successfully bound to the promoter to activate the expression of LacZ reporter gene. However, in the absence of TGA-element in Hap-2A-2 of TaPYL4-2A, TaARF4 could not bind the promoter to activate the LacZ reporter gene (Figure 7B). The results indicate that due to the difference of A/G at −1635 bp upstream of ATG translation starting point, TaARF4 could bind to the promoter region of Hap-2A-1, but not to the promoter region of Hap-2A-2.

We also performed a LUC experiment for further verification. Dual luciferase assays illustrated that LUC activity was higher than that of the negative control in the presence of TaARF4 effector and TaPYL4-2A promoter containing TGA-element (Figures 7C,D). On the contrary, LUC activity was similar to that of the negative control in the presence of TaARF4 effector and TaPYL4-2A promoter without the TGA-element. These results again indicated that TaARF4 functions as a transcription activator for Hap-2A-1, but not for Hap-2A-2.

According to the regional climate characteristics, topography, cultivation characteristics, sowing and maturity stage, Chinese wheat planting regions were divided into 10 major ecological areas (Zhang et al., 2002). From Chinese landraces to modern cultivars, the proportion of Hap-AB1 was increased significantly in regions I, II, VI, and IX, but not significantly changed in wheat regions III, IV, and V (Figures 8A,B). The mean PH in modern wheat varieties released in China from the 1940s to the 1990s has been reduced from 121.50 cm to 80.08 cm, and this change was associated with the increasing of proportion of haplotype Hap-AB1 from 16.70 to 72.41% (Figure 8C). To sum up, Hap-AB1 was selected by wheat breeders in the past decades, showing a considerable utilization potential in wheat genetic improvement.