A TP population containing 206 individuals was derived from a single F₁ seed from a cross between inbred line PH4CV as female parent and Zea mays ssp. mexicana (CIMMYT accession number 249743). In addition, an association panel including 351 elite inbred lines derived from a large association mapping population (Li et al., 2022) was used to conduct regional association mapping.

Genomic DNA from individuals in the TP population was extracted by the CTAB method. Genome sequencing libraries were sequenced using the Illumina Hiseq X platform, yielding a total of 1.18 Tb raw sequences (average depth 2.5×) with 150-bp paired-end reads. After checking and filtering sequence quality, the remaining sequences were mapped to the B73 reference genome (B73_V4). A total of 62,217,878 single nucleotide polymorphisms (SNPs) were identified in the population. To obtain high-quality SNPs, we removed low-quality SNPs with missing rate >60%, minor allele frequency (MAF) > 0.175 or MAF <0.075, and heterozygosity rate >5%. Finally, 138,208 high-quality SNPs were used to construct a recombination bin map.

The R package “binmapr” was used to fix genotypes using a windows size set to 15. Redundant markers were removed by the “bin” function in IciMapping software. The genetic distances between the final 1,951 bin markers were determined using a Kosambi method in IciMapping software with the “map” function. The linkage map and recombination bin map were drawn in R software with the “LinkageMapView” package and “plot” function, respectively.

To phenotype SRN, fifteen seeds of each line from the TP population and association panel were surface sterilized with 6% hypochlorite for 10 min, followed by three washes in distilled water, transferred to wet filter paper, rolled up and placed in plastic buckets containing distilled water. The buckets were incubated in darkness at 25°C. Three days later, the buckets were transferred to a new incubator at 23°C\8h darkness and 28°C\16h light. The distilled water in the buckets was replaced for every two days. After five days, the number of seminal roots for 10 healthy plants of each line was manually counted and average values of SRN were subsequently used in analysis.

QTL analysis of SRN based on linkage map and phenotypic values from the TP population was conducted using WinQTL Cartographer v2.5 software by composite interval mapping (CIM) (Wang et al., 2012), with a 1,000 permutations test at 95% confidence level to determine the optimal logarithm of odds (LOD) threshold values. A 1.5-LOD drop method was applied for defining the QTL confidence interval. Names were assigned to QTL following the nomenclature proposed by McCouch (1997), which combines the letter “q” for QTL, an abbreviation for the name of the trait, and a number for the chromosome.

Inbred lines IA2132 and PHW30 from the association panel were grown in a culture room. The number of seminal roots was determined from the seminal root primordia, which are initiated at the scutellar node of kernel about 25 days after pollination (Erdelska and Vidovencova, 1993). Samples from each inbred line for RNA-Seq were collected at the kernels scutellar node at the 20^th and 30^th days after pollination, respectively. Four biological replicates for each line were collected at each sampling time. All samples were immediately frozen in liquid nitrogen and stored at -80°C for RNA isolation. RNA extraction and library preparation for each sample were performed by Novogene Corporation. Sixteen libraries were sequenced using the Illumina HiSeq 2000 platform (Illumina, CA, USA), and 150 bp paired-end reads were generated. The raw reads were filtered to remove sequencing adapters as well as low quality reads (base number of Qphred ≤ 20 accounting for more than 50% of the read length). HISAT2 software was used to compare clean reads quickly and accurately with the reference genome to obtain the locations of reads on the AGv4 reference genome (Mortazavi et al., 2008). The FPKM (expected number of Fragments Per Kilobase of transcript per Million base pair sequences) was used to estimate gene expression levels (Trapnell et al., 2010). Differential expression analysis between the two inbred lines was performed via the pairwise comparison algorithm DESeq (Anders et al., 2013). DEGs were screened according to the following criteria: |log2 foldchange (FC)| ≥ 1, false discovery rate (FDR), and an adjusted P-value < 0.05 (Robinson et al., 2010). Gene Ontology (GO) enrichment analysis of the DEGs was conducted to identify the enriched biological functions between the two genotypes (Yu et al., 2012).

The genotype data from the association panel of 351 inbred lines in the qSRN7 interval were extracted from the resequencing data of 1,604 maize inbred lines reported by Li et al. (2022). After filtering using missing rates < 20%, MAF >0.05, a total of 33,164 high-quality SNPs were obtained. Principal component analysis (PCA) and kinship (K) matrix were calculated using 43,252 SNPs identified by the 50K SNP chip in the association panel. A mixed linear model (MLM) with PCA and K matrix was used to conduct regional association mapping in TASSEL 5.0 software (Bradbury et al., 2007). To determine the significant threshold for regional association results, we estimated the number of independent SNPs by pruning SNPs in the PLINK (window size 100, step size 50 and r² ≥ 0.2). After pruning, the number of independent SNPs in the qSRN7 interval was determined to be 6,070. We then selected 1.65×10^-4 (1/6070) as the threshold of association signals.

The TP population was genotyped using whole genome resequencing technology. A total of 138,208 high-quality SNPs were identified in the TP population and used to determine bin markers. Finally, a total of 1,951 recombination bins were obtained by using the “binmapr” package in R software ( Figure 1 ). The average physical interval of adjacent bins was 1,104.7 kb, with a maximum of 65,441.9 kb and minimum of 0.4 kb. Bin markers were mapped to maize B73 RefGen_V4 to assess the quality and accuracy of the map. The scatter plot of physical position of bin markers on all 10 chromosomes aligned well with the B73 reference genome, which indicated good collinearity between the maize B73 reference genome and the bin markers. The genetic length of the linkage map constructed using the 1,951 bin markers ( Supplementary Figure 1 ) was 1,142.1 cM with an average distance of 0.59 cM between adjacent markers. The number of markers per chromosome ranged from 102 (chromosome 9) to 329 (chromosome 4), with an average of about 195 bin markers per chromosome.

We observed that the SRN showed a large phenotypic difference between the PH4CV and teosinte parents ( Figure 2A ). The average SRN of the female parent PH4CV was 3.2, whereas the teosinte parent had no seminal roots. In the TP population developed from the cross between PH4CV and teosinte, SRN ranged from 0.3 to 4.5, with a coefficient of variation (C.V.) of 26.3% ( Figure 2B and Supplementary Table 1 ). Two QTLs were detected on chromosomes 1 and 7 ( Figure 2C ), explaining 6.18% and 6.25% of phenotypic variation, respectively. The known RTCS gene that has been cloned and regulates SRN, was located near the QTL peak on chromosome 1. The new QTL on chromosome 7 was named qSRN7, and a total of 345 candidate genes were located within the QTL interval.

We conducted a comparative transcriptome profiling of two inbred lines with extreme phenotypic differences in SRN. After filtration of low-quality sequences and adaptor sequences, a total of 670 million clean reads were obtained for the 16 RNA libraries. On average, 86.7% of reads were mapped uniquely to the B73 reference genome (AGPv4) ( Supplementary Table 2 ). Pearson correlation coefficients among the different biological replicates of the same genotype varied from 0.92 to 0.99 ( Supplementary Figure 2A ), suggesting high quality of the replicates. In the 16 samples, 59.8% of expressed genes were expressed at low levels (0 ≤ FPKM1< 1), 37.2% were expressed at medium levels (1 ≤ FPKM < 60), and 2.96% were expressed at high levels (FPKM ≥ 60) ( Supplementary Figure 2B ). To determine which expressed genes correlated with SRN development, DEGs were analyzed between line IA2132 with no SRN and PHW30 with multiple SRN at both sampling times ( Figure 3A ). Totals of 6,682 (3,114 down-regulated and 3,568 up-regulated) and 8,620 DEGs (3,397 down-regulated and 5,223 up-regulated) genes were identified at the 20^th and 30^th days after pollination, respectively ( Figure 3B ), with 4,241 DEGs shared in both samplings ( Figure 3C ). When compared with the qSRN7, 116 DEGs detected at least once were located within the qSRN7 interval ( Supplementary Table 3 ).

Among the DEGs, we successfully identified RTCS (Zm00001d027679) and BIGE1 (Zm00001d012883) at 30^th days after pollination, both genes regulate initiation of seminal roots in maize. We also identified promising DEGs associated with root development and growth. For example, ZmEXPB2 (Zm00001d029899) encoding an expansin-B4 protein, which increases primary root length (Zainab et al., 2021), was identified at both sampling times. Zm00001d029592 encoding a root hair defective (RHD3) protein, homologous to Arabidopsis RHD3 involved in regulating cell expansion and normal root hair development (Zheng and Chen, 2011), was also identified at .both times. GO analysis of shared DEGs identified at both sampling times revealed enrichment of genes in protein homodimerization activity, protein heterodimerization activity, heme binding, and tetrapyrrole binding ( Supplementary Figures 3A, B ). Previous studies documented that the interaction between RTCS and its paralogue RTCS-LIKE (RTCL) could form heterodimerization to regulate shoot-borne root initiation in maize (Majer et al., 2012), and RUM1 and its homolog RUM1-LIKE1 (RUL1) regulating seminal and lateral root formation can homo and heterodimerize in vivo (Zhang et al., 2016). These results suggested that those DEGs might be considered important candidate genes regulating seminal root formation in maize.

To identify candidate genes underlying qSRN7, we conducted regional association mapping in the association panel of 351 inbred lines. The panel showed large phenotypic differences in SRN ranging from 0 to 7 ( Figures 4A, B and Supplementary Table 4 ). Using the MLM model in TASSEL 5.0 software, we performed an association analysis of SRN within the qSRN7 interval. Five associated SNPs were detected at a threshold of −log10(P) > 3.78 ( Figure 4C ). We detected five candidate genes based on the physical location of those associated SNPs, including Zm00001d021861, Zm00001d021857, Zm00001d021572, Zm000001d021579 and Zm00001d021874.

To identify high-confidence candidate genes underlying qSRN7, we integrated DEGs obtained by RNA-Seq and regional association results. Three of the 5 candidate genes obtained by regional association analysis were differentially expressed between the inbred line with no SRN and the inbred line with multiple SRN; they were Zm00001d021572, Zm00001d021579 and Zm00001d021861. For Zm00001d021572 (encoding an anthocyanidin 3-O-glucosyltransferase), an associated SNP (S7_156440809) located within the UTR3 region, could form two haplotypes. Inbred lines with the GG haplotype had significantly more seminal roots than those with TT haplotype ( Figure 4D ). Morphologically, SRN of modern maize inbred lines differ substantially from their ancestors. To ascertain whether Zm00001d021572 underwent selection during maize domestication, we searched for this gene in selected regions by comparing maize and parviglumis using the XP-CLR method reported by Chen et al. (2022). We found that Zm00001d021572 was selected during maize domestication ( Figure 4E ). For Zm00001d021579 (encoding a DUF869 domain-containing family protein), an associated SNP (S7_156778214) located within the exon, could form two haplotypes. Inbred lines with TT haplotype had significantly more seminal roots than those with GG haplotype ( Figure 4F ). Zm00001d021861 encodes a vacuole membrane protein KMS1, which is involved in root development in Arabidopsis (Wang et al., 2011). An associated SNP (S7_ 164777805) located near the physical position of Zm00001d021861, could form two haplotypes. Inbred lines with CC haplotype had significantly more seminal roots than those with TT haplotype ( Figure 4G ). The expression data of these three candidate genes (Zm00001d021572, Zm00001d021579 and Zm00001d021861) in embryo, endosperm, seed, silks, tassel, cob, leaf, internode and root tissues were downloaded from MaizeGDB. The results showed that these genes had relatively high expression in the root tissues ( Supplementary Figure 4 ). These results suggested that all three candidate genes might play important roles in regulating SRN in maize.