Two corn varieties were used for hybridization in this study. The relatively short ear shank inbred line WL134 developed by our research team, featuring high comprehensive combining ability as the female parent. The long ear shank line L135 derived from a double haploid (DH) population constructed through the cross between WL134 and D7 inbred line with excellent resistance developed by our research team as the male parent. F₂ generations were harvested through self-pollination after F₁ planting. The F₂ mapping population and parents were planted during the summer of 2024 at Dishang Experimental Station (114°43′7.928″ E, 37°56′25.800″ N) of the Institute of Cereal and Oil Crops, Hebei Academy of Agriculture and Forestry Sciences. Parents were planted 2 rows and F₂ population planted 25 rows. The plants were grown with a row length of 6 m, row spacing of 0.5 m, under standardized production conditions for planting and management. At maturity, the ear shank length (actual length from the bottom of the ear to the stem attachment point) was measured with band tapes repeating three times. Due to the curved structure of the ear shank, the tape was kept closely adhered to the ear shank during measurement to minimize errors.

Through field cultivation and measurement of ear shank length, 30 plants with extremely short ear shanks were selected from the F₂ population to form the short ear shank bulk (S-pool), and 30 plants with extremely long ear shanks were selected to form the long ear shank bulk (L-pool). The corresponding individual plants were numbered, and leaf samples were collected from each for preliminary mapping of the maize ear shank length genes.

BSA-seq was used to identify the genes regulating ear shank length in the F₂ population. We selected 30 plants that were extremely long and short ear shank length to create an extreme population. Leaf samples were collected from three individuals of each parental line, as well as from the 30 plants in each of the extreme ear shank bulks (S-pool and L-pool). These samples were submitted to SMART GENOMICC (Qingdao, China) for BSA-seq analysis. The sequencing depth for the parental lines was set at 20×, while that for the offspring pools was set at 30×.DNAsecure Plant Kit (TIANGEN) was used for DNA extraction from plant tissues. The DNA quality was checked using NanoDrop 2000 spectrophotometer (Thermo Fischer Scientific), Agarose gel electrophoresis and Qubit fluorometer (Invitrogen). Qualified DNA samples were randomly broken into 350bp fragments by a Covaris crusher. Sequencing libraries were constructed following the segmentation of DNA through end repair, addition of adenine to the 3’ ends, adapter ligation and PCR amplification. Once the library passed quality checks, sequencing was performed using the DNBSEQ-T7 platform. The default parameters of fastp (Chen et al., 2018) were used for quality control, unqualified reads were filtered, and the clean reads obtained were used for subsequent analysis. Clean reads were aligned to the reference genome sequences of the B73 genome (https://download.maizegdb.org/Zm-B73-REFERENCE-NAM-5.0/) using Sentieon software (Pei et al., 2021). The mapping results were sorted and deduplicated using Samtools (Li et al., 2009). SNPs and InDels were detected and annotated using Sentieon (Pei et al., 2021) and ANNOVAR software (Wang et al., 2010). The SNP/InDels were selected with read depth ≥ 4, RMS Mapping Quality ≥ 40 and Genotype Qualit ≥ 5. The association analysis was conducted with ΔSNP/InDel index (Fekih et al., 2013) and G′-value (Magwene et al., 2011). The overlapped regions based on the above two methods were considered candidate regions for ear shank length.

The collected samples from three replicated plants were stored in refrigerator at -80°C. The ear shank at silking stage of WL134 and L135 were used for transcriptome sequencing. DNA library generation and RNA-seq high throughput sequencing were performed by SMART GENOMICC (Qingdao, China). The total RNA was extracted from the samples using standard extraction methods. RNA quality control was assessed with the Agilent 2100 bioanalyzer. The mRNA with a polyadenylic acid tail was enriched by connecting oligothymidine magnetic beads, and then the obtained mRNA was randomly interrupted with divalent cations in NEB fragmentation buffer (Parkhomchuk et al., 2009). mRNA quality, including the mRNA concentration and fragment size, was tested by using Qubit2.0 and Agilent 2100. A total of 6 qualified libraries were sequenced on the Illumina Novaseq platform HiSeqTM 2500 (Illumina, San Diego, CA, USA), and 150 bp paired-end reads were generated.

Raw sequence reads were processed using fastp (Chen et al., 2018) for quality control to filter the unqualified reads with default parameters. The clean reads were mapped to a B73 reference of maize genome (https://download.maizegdb.org/Zm-B73-REFERENCE-NAM-5.0/) using hisat2 (Mortazavi et al., 2008). The number of fragments per kilobase of transcript per million mapped reads (FPKM) value for each gene was calculated using Featurecounts (Liao et al., 2014). Differential expression analysis between the two comparative combinations was performed using DESeq2 software (Love et al., 2014). The method of Benjamini and Hochberg was used to adjust the resulting P-values to control for false discovery rates. Genes with adjusted P-values <0.05 were classified by DESeq2 as differentially expressed. Gene Ontology (GO) divides the functions of genes into three parts: cellular component (CC), molecular function (MF) and biological process (BP). The ClusterProfiler software was used to perform GO functional enrichment analysis and Kyoto Encyclopedia of Genes and Geno (KEGG) pathway enrichment analysis on the differential gene sets (Wu et al., 2021). Padj less than 0.05 was used as the threshold of significant enrichment.

Total RNA was extracted from the ear shanks at silking stage using an RNA pure plant kit (TIANGEN, Beijing, China). The first strand of cDNA was synthesized according to the instructions for the HiScript^®III RT SuperMix for qPCR (Vazyme, Nanjing, China). Subsequently, qRT‐PCR was performed using the ChamQ SYBR Color qPCR Master Mix (Vazyme, Nanjing, China). qRT-PCR was performed on CFX384 Real-Time System (BIO-RAD, United States). The maize Tubulin gene was used as the internal control. The 2^–ΔΔCT (Livak and Schmittgen, 2001) quantitative analysis method was used to calculate the relative expression level. The primers used in this study are listed in Supplementary Table S1.

Observing the phenotypic traits related to ear shank length in the two parental lines WL134 and L135, the average ear shank length of WL134 was 7.97 cm, while that of L135 was 21.83 cm. WL134 represents a relatively short ear shank material, with its ear shank length showing a significant difference compared to L135 (Figures 1A, B). Additionally, ear shank related traits including plant height, ear height and internodes number were examined in both parents. Significant differences were observed across these traits between WL134 and L135 (Figures 1C-E).

To map the genes regulating ear shank length, we constructed an F₂ population with 308 lines derived from a cross between the long ear shank length inbred line L135 and the relatively short ear shank length inbred line WL134. This F₂ population was used to identify phenotypic characterization of ear shank length. The F₂ plants were phenotyped, 30 lines exhibiting long ear shank length (such as L135 and L1-30) and 30 individuals exhibiting short ear shank length (such as WL134 and S1-30) were selected based on ear shank length measurements to construct the corresponding extreme bulks (Supplementary Figure S1). The survey and analysis of ear shank length revealed a significant difference between the two extreme pools, which was suitable for follow-up research (Figure 1F). So, these selected lines were used to construct the extreme long ear shank length bulk and the extreme short ear shank length bulk for BSA-seq analysis.

Sequencing was performed on the parental materials as along with the extreme long ear shank bulk and short ear shank bulk derived from the F₂ population (Table 1). A total of 31.59 million clean reads were generated for the relatively short ear shank parent WL134, 28.73 million clean reads were generated for the long ear shank parent L135. 40.47 million clean reads were generated for long ear shank bulk, and 45.48 million clean reads were generated for short ear shank bulk. The sequencing depths for WL134, L135, long ear shank bulk and short ear shank bulk were 22.63×, 20.47×, 25.51× and 30.38×, respectively. The properly paired ratios were more than 87.51% (97.40% for WL134, 96.01% for L135, 87.51% for L-pool, 93.13% for S-pool). All samples exhibited high data alignment rates, providing a reliable foundation for subsequent SNP detection.

A total of 1,503,803 SNPs and 197,389 Indels were finally detected by comparison and filtering between the two bulks (Supplementary Figure S2). Association analysis between the ear shank length trait and polymorphic markers were performed using the Δ(SNP/InDel-index) and G′-value methods. After overlapping the results obtained based on the 99% significance level (top 1%, red line) thresholds of the Δ(SNP/InDel-index) confidence intervals (Figure 2A) and the G′-value thresholds (Figure 2B), fourteen QTLs were identified (Supplementary Table S2). These loci were located on Chr1, Chr2, Chr4, Chr5, Chr6, Chr7, Chr8 and Chr9. Among them, qESL1, qESL4, qESL5, qESL7 and qESL10, five loci co-localized with previously reported QTLs, whereas the remaining nine represent novel genomic regions associated with ear shank length.

To investigate the DEGs between the long ear shank line L135 and the relatively short ear shank line WL134, we performed transcriptome sequencing analysis on different ear shank samples at the silking stage. The RNA sequencing of six cDNA libraries (L135-1, L135-2, L135-3, WL134_1, WL134_2, WL134_3) were generated after filtering a total of 44.87 Gb clean bases (Table 2). The average percentages of Q20 and Q30 were 98.0% and 94.5%, respectively. 91.96% to 92.66% of the clean reads were successfully mapped to the reference genome B73_V5 using the HISAT2 software. Pearson correlation analysis revealed that the correlation between replicates was stronger than that different samples, further demonstrating the accuracy and reproducibility of the sequencing results (Supplementary Figure S3). A total of 3,460 DEGs were identified between the long ear shank and the short materials. Compared to the short ear shank material WL134, 1,802 genes were upregulated and 1,658 genes were downregulated in the long ear shank material L135 (Figure 3A).

KEGG enrichment analysis revealed that all the DEGs were enriched in 113 metabolic pathways. The top 20 KEGG pathways were based on enrichment factor (Figure 3B). The DEGs were primarily involved in pathways such as phenylpropanoid biosynthesis, zeatin biosynthesis, glyoxylate and dicarboxylate metabolism, brassinosteroid biosynthesis and pyruvate metabolism, etc. The pathways including flavonoid biosynthesis, phenylalanine metabolism, tropane, piperidine and pyridine alkaloid biosynthesis and tyrosine metabolism were exclusively enriched among upregulated genes, while the remaining pathways were contained both upregulated and downregulated genes.

GO analysis was conducted to explore the BP, CC and MF associated with ear shank length variation between WL134 and L135. The top 10 GO terms were selected based on the number of enriched genes (Figure 3C). Common GO terms included movement of cell or subcellular component, microtubule-based movement, response to oxidative stress, external encapsulating structure organization and cell wall organization were enriched in BP. GO terms included extracellular region, cell wall, external encapsulating structure, apoplast and photosystem II were significantly enriched in CC. GO terms included heme binding, tetrapyrrole binding, incorporation or reduction of molecular oxygen, iron ion binding and protein dimerization activity were notable enriched in MF.

To further explore candidate genes associated with maize ear shank length, the results obtained by BSA-seq and RNA-seq were combined to analyze. Based on polymorphisms identified between the parents, we detected 334 genes containing nonsynonymous mutations, synonymous mutations and frameshift mutations in their open reading frames (ORFs) regions by BSA-seq. We compared the 334 genes obtained by BSA-seq with the 3460 DEGs obtained by RNA-seq, only nineteen common genes were selected. Through annotation analysis of SNPs within genes (Supplementary Table S3), it was found that the Zm00001eb166380 gene had a frameshift insertion, while the Zm00001eb023420, Zm00001eb166550 and Zm00001eb353540 genes had frameshift deletions. Stopgain SNVs were present in the Zm00001eb353540 and Zm00001eb382210 genes. Synonymous SNVs, nonsynonymous SNVs, nonframeshift deletions or nonframeshift insertions were present in the exons of the other genes. The RNA-seq results revealed differential expression of these 19 genes between the two parental lines (Figure 4). Among these, five genes (Zm00001eb023280, Zm00001eb282410, Zm00001eb282430, Zm00001eb353170 and Zm00001eb353540) were enriched in BP ontologies. One gene (Zm00001eb023280) was enriched in CC ontologies, and six genes (Zm00001eb023280, Zm00001eb166400, Zm00001eb282410, Zm00001eb282430, Zm00001eb316880 and Zm00001eb353540) were enriched in MF ontologies by GO function (Supplementary Table S3).

To further confirm whether the candidate genes exhibit differential expression levels between the two parental lines and to validate the accuracy of the transcriptomic differential expression data, associated DEGs identified through the integrated BSA-seq and RNA-seq analysis were selected for validation using quantitative real-time PCR (qRT-PCR). These results indicated that the expression trends of these genes aligned consistently with the RNA-seq results (Figure 5). Among them, six genes including Zm00001eb023400, Zm00001eb023420, Zm00001eb166550, Zm00001eb316880, Zm00001eb282410 and Zm00001eb282430 were significantly upregulated in the long ear shank lines. Conversely, seven genes including Zm00001eb117860, Zm00001eb117870, Zm00001eb243020, Zm00001eb050490, Zm00001eb166380, Zm00001eb316700 and Zm00001eb353170 were significantly downregulated (Table 3). Consequently, the integration of BSA-seq and RNA-seq analyses identified 13 genes as high-confidence candidates within the targeted genomic regions.

Zm00001eb023400 and Zm00001eb023420 on chromosome 1 exhibited upregulation, whereas Zm00001eb050490 showed downregulation in comparison to the short ear shank length line WL134 (Table 3). All three genes encoded lytic transglycolase situated in the membrane, primarily functioning as components of the cell wall or cellular component curator. On chromosome 2, Zm00001eb117860 and Zm00001eb117870 were downregulated. The two genes belonged to ATPase family and regulated the binding of sequence specific DNA. Based on the annotation information, Zm00001eb166380 and Zm00001eb166550 on chromosome 4 encoded ring domain ligase 2 and glycine rich domain containing protein respectively. Zm00001eb243020 on chromosome 5 encoded protein smax1-like 3. On chromosome 6, Zm00001eb282410 and Zm00001eb282430 encoded peroxidase, all of which were enriched in phenylpropanoid biosynthesis pathway. Zm00001eb316700 and Zm00001eb316880 on chromosome 7 encode IRK-interacting protein and wall associated receptor kinase respectively. Zm00001eb353170 on chromosome 8 encoded non-specific serine/threonine protein kinase.