The 254 RIL-F_2:8 population were obtained by single seed decent (SSD) from the F₂ population of a cross between parental inbred lines B73 and Abe2. The representative inbred line B73 exhibited a leaf scorching (thermosensitive) phenotype when subjected to high temperature stress above 35°C during the flowering stage in a field environment. And the native waxy maize inbred line Abe2 in northwestern China exhibited a leaf scorch-free (thermotolerant) phenotype under the same conditions. The 254 RIL-F_2:8 population and their parents used to construct the high-density genetic map were planted in Ledong County, Hainan, China (18°45N, 109°10E) in 2016.

Research by Frey et al. showed that leaf scorching trait was one of the important thermosensitive phenotypes under high temperature stress at flowering in maize (Frey et al., 2016). The leaf scorching phenotype in this study was divided into three categories: Leaf Scorching damage (LS), Leaf Scorching Degree (LSD), and Leaf Scorching Ratio (LSR). Phenotypic data of the three types for the LS trait were collected by visual method: extreme leaf scorching damage (phenotypes of leaf scorching are consistent with B73), leaf scorch-free damage (phenotypes of leaf scorching are consistent with Abe2.), and intermediate type. LSD represented the proportion of the leaf scorching damage area in the whole leaves, indicating the degree of leaf scorching damage. It is indicated by the Roman number 1 (no leaf scorching damage) to 9 (extreme leaf scorching damage), with a total of 9 indication levels. LSR referred to the ratio of the number of scorching leaves to the total number of leaves, expressed as 0 to 100%. The broad-sense heritability of leaf scorching traits was calculated according to the following formula:

where σ_g ² represented the variance of genetic effects, and σ_e ² represents the variance of environmental effects (Yang et al., 2019). A basic statistical analysis was implemented by the SPSS16.0 software with default parameters (SPSS Inc., Chicago, IL, USA) (Zhu et al., 2019).

In view of the continuous high temperature weather conditions in Hefei over the past years, we planted the 254 RIL-F_2:8 population and their parents for three biological replicates per year in Dayang experimental farm of Anhui Agricultural University (31°49N, 117 °13E) with interval contrast design (ICD) in early June 2017 and 2018 (Zeng et al., 2020). During the maize jointing to flowering stage, they were subjected to continuous high temperature stress, as shown in Supplementary Figure 1 . In fact, 24 and 29 days of high temperatures above 35°C in 2017 and 2018, respectively.

When the 254 RIL-F_2:8 population and parental inbred lines grew to the V6 period (6 visible leaves), fresh green leaves were removed and stored in the dry ice. Extraction of total genomic DNA for experimental samples using a modified Cetyltrimethyl ammonium bromide (CTAB) method (Murray and Thompson, 1980). An improved SLAF-seq strategy was utilized in our experiment. In this experiment, Oryza sativa L Geng/japonica was used as the control, and the evaluation of the control data monitored whether the experimental process was normal and the effectiveness of the digestion protocol was determined. First, the pre-designed scheme of SLAF is selected using the training data. B73_RefGen_V4 reference genome of maize was used to simulate the number of markers produced by different enzymes, designing marker identification experiments. Next, SLAF library construction was conducted according to a pre-designed scheme. For the 254 RIL population, two enzymes (HaeIII and Hpy166II, New England Biolabs, NEB, USA) were used to digest the genomic DNA at 37°C. A single nucleotide (A) overhang was added subsequently to the digested fragments using Klenow Fragment (3´→ 5´) and dATP (New England Biolabs, NEB, USA). Duplex tag-labeled sequencing adapters (PAGE-purified, Life Technologies, USA) were then ligated to the A-tailed fragments using T4 DNA ligase. Polymerase chain reaction (PCR) was performed using diluted restriction-ligation DNA samples, dNTP, Q5^® High-Fidelity DNA Polymerase and PCR primers (Forward primer: 5’-AATGATACGGCGACCACCGA-3’, reverse primer: 5’-CAAGCAGAAGACGGCATACG-3’) (PAGE-purified, Life Technologies). PCR products were then purified using Agencourt AMPure XP beads (Beckman Coulter, High Wycombe, UK) and pooled. Pooled samples were separated by 2% agarose gel electrophoresis. Fragments ranging from 414 to 464 base pairs (with indexes and adaptors) in size were excised and purified using a QIAquick gel extraction kit (Qiagen, Hilden, Germany). Gel-purified products were then diluted. And pair-end sequencing was performed on an Illumina platform system (Illumina, Inc; San Diego, CA, USA) according to the manufacturer’s recommendations.

The SLAF marker identification and genotyping were performed according to procedures described by Sun et al. (Sun et al., 2013). Briefly, low-quality reads were filtered out and then raw reads were sorted to each subsequence according to duplex barcode sequences. After the barcodes and the terminal 5-bp positions were trimmed from each high-quality reads, clean reads from the same sample were mapped onto the maize genome sequence using SOAP software (Li et al., 2008). The reference genome information is based on version B73_RefGen_v4. Sequences mapping to the same position were defined as one SLAF locus (Zhang et al., 2015). Single nucleotide polymorphism (SNP) loci of each SLAF locus were then detected between parents, and SLAFs with more than 3 SNPs were filtered out firstly. In order to obtain high-quality SLAF markers for genetic map construction, first, the average sequence depth should be >21X for parents, while for each offspring the reads with sequence depth >8X were used to define alleles. Second, markers with a data missing rate of more than 50% were filtered. Third, Chi-square test was used to detect segregation distortion. In the process of map construction, markers with significant partial separation (P<0.05) were treated as an auxiliary marker. SLAFs with two to four alleles were identified as polymorphic and considered potential markers. All polymorphism SLAFs loci were genotyped with consistency in the parental and offspring SNP loci. The marker code of the polymorphic SLAFs were analysed according to the RIL population type, which consisted of one segregation types (aa×bb).

SLAF loci were partitioned primarily into linkage groups (LGs) based on their locations on B73_RefGen_V4 reference genome of maize. Next, the modified logarithm of odds (MLOD) scores between SLAF markers were calculated to further confirm the robustness of markers for each LGs. Markers with MLOD scores< 5 were filtered prior to ordering. To ensure efficient construction of the high-density and high-quality map, a newly developed HighMap strategy was utilized to order the SLAF markers and correct genotyping errors within LGs (Liu et al., 2014). Firstly, recombinant frequencies and LOD scores were calculated by two-point analysis, which were applied to infer linkage phases. Then, enhanced Gibbs sampling, spatial sampling and simulated annealing algorithms were combined to conduct an iterative process of marker ordering (Jansen et al., 2001). Briefly, in the first stage of the ordering procedure, SLAF markers were selected using spatial sampling. One marker was taken randomly in a priority order of test cross, and markers with a recombination frequency smaller than a given sampling are excluded from the marker set. Subsequently, simulated annealing was applied to searching for the best map order. Summation of adjacent recombination fractions was calculated as illustrated by (Liu et al., 2014). The annealing system continued until, in a number of successive steps, the newly generated map order is rejected. Blocked Gibbs sampling was employed to estimate multipoint recombination frequencies of the parents after the optimal map order of sample markers were obtained. The updated recombination frequencies were used to integrate the two parental maps, which optimize the map order in the next cycle of simulated annealing. Once a stable map order was obtained after 3-4 cycles, we turned to the next map construction round. A subset of currently unmapped markers was selected and added to the previous sample with decreased sample radius. The mapping algorithm repeats until all the markers were mapped appropriately. The error correction strategy of SMOOTH was then conducted according to parental contribution of genotypes (Van Os et al., 2005), and a k-nearest neighbor algorithm was applied to impute missing genotypes (Huang et al., 2021). Skewed markers were then added into this map by applying a multipoint method of maximum likelihood. Map distances were estimated using the Kosambi mapping function (Kosambi, 2016).

The QTL mapping of thermotolerant traits was performed by R/qtl software (Broman et al., 2003) for composite interval mapping (CIM) analysis, and the logarithm of odds (LOD) significance threshold levels was determined by 1000-permutation test (P<0.05). The confidence interval for each QTL was defined using a 2-LOD support interval (Li et al., 2016). The mapping interval of each QTL was determined by the the peak of the LOD and its surrounding value (≥ 2). The software used for drawing the map were origin 7.0 and HighMap. We estimated additive effects and the phenotypic variance explained by individual QTL by the coefficient of determination (R²). Positive additive effects indicated favorable alleles derived from Abe2, while negative additive effects indicated favorable alleles from B73.

After the major QTLs for LS trait were identified, the genetic effect analysis of each QTL locus was carried out to determine the QTL locus with the greatest effect. Combined with previously published QTL-seq results to identify final candidate genes controlling thermotolerant traits (Zeng et al., 2020). The qPCR primers of candidate genes were designed to carry out relative quantitative analysis, with ubiquitin as an internal reference gene. ( Supplementary Table 1 ). The parental lines Abe2 and B73 were germinated and grown in a growth chamber at 28°C until the stage of visible 4-leaf growth. Subsequently, a high-temperature stress experiment was conducted in a plant light incubator at a temperature of 42°C and a humidity of 65%. 2×Spark Taq PCR Master Mix (with dye) (Shandong Sparkjade Biotechnology Co., Ltd.) is used for PCR amplification. TRIzol kit (TIANGEN, W9330) was used to extract RNA from fresh leaf, and SPARKscript II RT kit (With gDNA Eraser) (Shandong Sparkjade Biotechnology Co., Ltd.) was used to reverse message RNA into cDNA, which was diluted by five times and amplified by 2×SYBR Green qPCR Mix (With ROX) (SparkJade, AH0104-B).

To develop high-density polymorphic molecular markers, we performed high-throughput SLAF sequencing of the 254 RIL-F_2:8 population. This resulted in a total of 1,366.12 million reads and 273.03 Gb raw data, with an average Q30 percentage of 93.03% and an average GC percentage of 42.79% ( Supplementary Figure 1A, B and Table 1 ). For the paternal inbred line (Abe2), a total of 684,620 SLAF markers were developed with an average sequencing depth of 18.48×. The maternal inbred line (B73) produced 955,168 SLAF markers with an average sequencing depth of 23.95×. In addition, 420,299 SLAF markers were developed for the RIL-F_2:8 population, with an average sequencing depth of 8.05× ( Supplementary Table 2 ).

The SLAF markers were aligned to the maize reference genome using BWA software (Van Os et al., 2005), and the count of SLAF markers and polymorphic SLAF markers on each linkage group was determined. The length of SLAF markers ranged from 414 bp to 464 bp ( Supplementary Figure 1C ). Polymorphism analysis was conducted based on the variation in allele number and marker sequence, resulting in three types of SLAF markers: Polymorphic, Non-Polymorphic, and Repetitive SLAF ( Figure 1 ; Supplementary Table 3 , 4 ). A total of 589,770 SLAF markers developed for all samples were used for SLAF polymorphism marker analysis ( Supplementary Table 3 ). Among them, 108,709 were identified as polymorphic SLAF markers, with a polymorphism ratio of 18.43% ( Supplementary Table 3 ). The distribution of SLAF markers and polymorphic SLAF markers across each linkage group was visualized ( Figures 1A, B ; Supplementary Table 4 ), indicating an even distribution of SLAF markers on each linkage group.

To facilitate further genetic analysis, the polymorphic markers were encoded following the general 2-allelic coding rules of genetics ( Supplementary Table 5 ). After filtering out the SLAF markers with parental information, a total of 108,709 polymorphic SLAF markers were genotyped and classified into eight segregation patterns ( Figure 1C ). Additionally, we calculated the recombination values between all polymorphic SLAF markers and observed fewer recombination hotspots on each chromosome ( Figure 1D ). Based on the genetic characteristics of the RIL-F_2:8 population, the aa x bb type polymorphic markers were selected as valid markers. Consequently, 48,397 markers belonged to the valid marker category of the aa x bb type, resulting in an effective polymorphism rate of 8.21% ( Supplementary Table 3 ).

To ensure the quality of the genetic map, the polymorphic SLAF markers underwent several filtering steps. First, markers were filtered based on parental sequencing depth, ensuring it was above 6X. Then, markers with less than 5 SNPs were excluded. Completeness filtering was applied, requiring the genotype to cover at least 50% of all progeny individuals. Additionally, markers showing partial separation were filtered out. After this process, 10,724 polymorphic SLAF markers remained for genetic map construction. Subsequently, markers with MLOD values lower than 5 when compared with other SLAF markers were further filtered (Huang et al., 2011), resulting in a final count of 10,112 polymorphic markers used for constructing the genetic map ( Table 2 ).

Each chromosome served as a linkage group, and HighMap software (Li et al., 2018) was utilized to analyze the linear array of SLAF markers. Genetic distances between adjacent markers were estimated, culminating in the construction of a genetic map with a total map distance of 1,475.88 cM ( Table 2 ; Figure 2 ). The high-density genetic map had an average genetic distance of 0.16 cM, with 98.77% of the gaps being less than 5 cM and distributed nearly evenly across each linkage group ( Table 2 ). Among the 10 chromosomes, chromosome 9 exhibited the largest gap of 10.9 cM due to an insufficient number of polymorphic markers, resulting in an average genetic distance of 0.73 cM, the highest among all linkage group maps. Conversely, chromosome 5 boasted the highest number of markers, with 2,121 polymorphic SLAF markers and the smallest average genetic distance of 0.09 cM among all linkage group maps ( Table 2 ).

In order to assess the quality of the genetic map, a comprehensive evaluation of various aspects was conducted to ensure its accuracy and reliability. These aspects include analyzing the collinearity between the map position of SLAF markers and their physical location on the genome ( Supplementary Table 6 ), examining the sequencing depth of the markers ( Supplementary Table 7 ), analyzing SNP markers within each linkage group ( Table 3 ), assessing segregation distortion ( Supplementary Table 9 ), evaluating the integrity of markers across all individuals ( Supplementary Figure 3 ), conducting monomer analysis for individual genotypes ( Supplementary Figure 4 ), and studying the recombination relationship between SLAF markers and adjacent ones ( Supplementary Figure 4 ). Each aspect provided valuable insights into the quality and reliability of the genetic map.

The quality of the genetic map has a direct impact on subsequent QTL mapping and genetic analysis. Upon evaluating the seven key aspects mentioned above, it is evident that this experiment yielded a high-quality, high-density genetic map. The map order of SLAF markers aligns consistently with the physical genome location ( Supplementary Table 6 ), with each molecular marker exhibiting an average sequencing depth exceeding 10X ( Supplementary Table 7 ). Furthermore, the genotype integrity of each individual reached 99.83% ( Supplementary Figure 3 ), and the allelic origins for larger genetic segments in each individual demonstrated high consistency ( Supplementary Figure 4 ). In summary, we utilized SLAF-seq technology to develop 10,112 polymorphic molecular markers, resulting in the construction of a high-quality, high-density genetic map.

Frey et al. proposed that leaf scorching is a crucial thermosensitive phenotype during the reproductive stage of maize under high temperature stress (Frey et al., 2016). A survey of daily temperatures throughout the maize growth period revealed that the parental inbred lines and RIL-F_2:8 population encountered temperatures exceeding 35°C for up to 17 days in 2017 and 19 days in 2018 ( Supplementary Figure 1 ). Following prolonged exposure to high temperatures, the parent B73 displayed noticeable leaf scorching damage, whereas the parent Abe2 remained unscathed by leaf scorching (Zeng et al., 2020). The RIL-F_2:8 population exhibits a remarkably diverse range of phenotypic variations, characterized by various levels of thermosensitive phenotypes under prolonged exposure to high temperature stress. ( Supplementary Figure 5 , 6 ; Supplementary Table 8 ). Correlation coefficients between the three indicators of the leaf scorching trait were calculated as 0.146 (between LS and LSR), 0.144 (between LS and LSD), and 0.875 (between LSR and LSD) ( Supplementary Figure 6 ). Notably, a highly significant positive correlation between LSR and LSD was observed, suggesting these two indicators could serve as reliable references. The broad-sense heritability values for the three indicators were 0.73 (LS), 0.38 (LSD), and 0.60 (LSR), respectively ( Supplementary Table 10 ). Despite the strong correlation between LSR and LSD, their heritability differed significantly, with LS exhibiting the highest heritability and LSD the lowest. Overall, this research sheds light on the complex relationship between high temperature stress and leaf scorching in maize, providing valuable information for understanding and potentially improving the plant’s resilience to such environmental challenges.

Based on the high-density genetic map and phenotypic characterization, we performed QTL analysis of thermotolerance traits using composite interval mapping (CIM) method of R/qtl software. A total of 16 QTLs were identified ( Table 4 , Supplementary Table 11 ; Figure 3 ), including 4 QTLs belonging to LS, and 6 QTLs each for LSD and LSR, respectively. The four LS QTLs were distributed on chromosomes 1 (qLS1), 2 (qLS2.1, qLS2.2) and 3 (qLS3), explaining a phenotypic variation range between 2.83% and 6.58% ( Table 4 ). The six LSD QTLs were located on chromosome 1, 2, 5, 6, 8 and 9 respectively, explaining a phenotypic variation range between 2.13% and 9.11%. The genetic architecture of LSD featured a large-effect QTL along with many small-effect QTLs ( Table 4 ). The six LSR QTLs were located on chromosomes 1, 4, 5, 6, and 8, respectively, collectively explaining 25.98% of the phenotypic variation. Further genetic architecture information for thermotolerance traits could be found in ( Table 4 , Figure 3 ; Supplementary Figure 7 ). Among the 16 QTLs, the favorable allele loci of 11 QTLs were from the parental inbred line Abe2, while the favorable allele loci of the remaining QTLs were from the maternal inbred line B73 ( Supplementary Figure 7 ). The localization results and chromosomal distribution of the 16 QTL loci associated with thermosensitive phenotype traits are presented in Supplementary Figure 8 . Chromosome 1, 2, and 8 each harbor three QTL loci, while chromosome 5 and 6 contain two QTL loci each. Lastly, chromosome 3 and 4 each have one QTL locus. Most of the QTL loci span large genomic regions.

Among the 16 identified QTLs associated with thermotolerance traits, the qLS1 allele displayed a notable difference of 8.13E-05 ( Figure 4A ), pointing to its significance as a key candidate locus. Additionally, we conducted an analysis of the genetic effects related to four QTL loci governing the LS trait, with the outcomes depicted in Figure 4B . Notably, when the allele of qLS2.2 was denoted as a1, and the alleles of the other three loci were a1 or b1, the thermotolerance index of the RIL-F_2:8 individuals generally exhibited lower levels, highlighting the pivotal role of the qLS2.2 locus in controlling the LS trait. Furthermore, the majority of RIL-F_2:8 individuals displayed leaf scorching phenotypes when possessing the a1 allele of the qLS2.2 locus. Conversely, when the allele of qLS2.2 was b1, and alleles of qLS2.1 and qLS3 were a2a3, the RIL-F_2:8 individuals demonstrated higher thermotolerance indices, particularly when the allele of qLS1 was a4, signifying its crucial role in regulating the LS trait ( Figure 4B ). Notably, the allele of qLS1 was a4, the RIL-F_2:8 individuals could exhibit strong thermotolerance. Furthermore, in conjunction with preceding QTL-seq results from our laboratory (Zeng et al., 2020), a co-localization on chromosome 1 was discovered, encompassing only six genes governing the thermotolerance trait ( Figure 4C ).

The annotation details for the six candidate genes are provided in Table 5 . Gene 1 is identified as a 2-oxoglutarate (2OG) and Fe (II)-dependent gene, potentially implicated in redox processes. Gene 2 is characterized as an N-acetyl-gamma-glutamyl-phosphate reductase, also likely involved in redox processes. Gene 3 belongs to the Cyclin class and may play a role in regulating the cell cycle. Gene 4 is classified as an ARM repeat superfamily protein with no associated functional annotation. Gene 5 is described as a Calcium-transporting ATPase 2 plasma membrane-type, possibly involved in ion transport processes. Lastly, Gene 6 is identified as a LOB domain-containing protein 28, potentially involved in ATP hydrolysis-coupled cation transmembrane transport processes.

Subsequently, the coding sequence variants of the six candidate genes were analyzed, revealing that only four genes exhibited variations in the coding region ( Supplementary Table 12 ). Notably, gene 2 harbored a sole non-synonymous variant SNP. Analysis of the tissue expression patterns of the six candidate genes indicated robust expression of gene 2, particularly in leaves and internodes, while the remaining five genes displayed markedly lower expression levels ( Figure 4D ). Furthermore, assessment of the relative expression of the six candidate genes under high-temperature stress unveiled significant differences in the expression of gene 2 between the parents (P<0.05) ( Figure 4E ), with gene 5 exhibiting low expression without significant variance ( Figure 4F ). Conversely, the expression of the remaining four genes was undetectable due to extremely low expression levels. In summary, based on significant differences in both expression levels and coding region sequences, gene 2 emerges as the prime candidate gene influencing the thermotolerance trait.

We further analyzed the variations in the promoter region and found a total of 1681 SNP variations within the QTL mapping interval, with 22 SNPs located in the coding sequence (CDS) region, most of which were promoter region variations ( Supplementary Table 13 ). Specifically, the promoter region of Zm00001d033328 exhibited 47 SNP variations along with 41 cis-regulatory elements ( Supplementary Table 14 ). These elements play a crucial role in various biological processes such as MeJA-responsiveness, light responsiveness, regulation of flavonoid biosynthetic genes, auxin-responsive elements, gibberellin-responsive elements, drought-inducibility, abscisic acid responsiveness, and MeJA-responsiveness.