In this study, 380 maize grain samples representing 190 different maize inbred lines were selected as experimental materials. These samples encompass a wide range of genetic backgrounds and originate from China, Thailand, Mexico, and Peru, among other locations. For each maize inbred line, two samples were selected, each consisting of 5 kernels for grinding and subsequent folate determination using the HPLC method. Based on the median folate content, the maize inbred lines were categorized into high-folate and low-folate groups. The sample details have been provided in Supplementary Table S1 .

The total RNA was extracted using the RNAprep Pure Plant Kit (Tiangen, Beijing, China) according the instructions provided by the manufacturer. The constructed libraries were sequenced on the Illumina Novaseq 6000 sequencing platform. The raw reads were further processed with a bioinformatic pipeline tool, BMKCloud (www.biocloud.net) online platform. To ensure the accuracy, reads with more than 10% N bases and low-quality reads with Q ≤ 10 and over 50% bases were excluded. After filtering low-quality reads, these clean reads were then mapped to the reference genome sequence (B73 RefGen_v3) using Hisat2. Only reads with a perfect match or one mismatch were further analyzed and annotated based on the reference genome.

Differential expression analysis of two groups was performed using the DESeq2. Gene expression levels were calculated using fragments per kilo base of transcript per million mapped reads (FPKM). The resulting P values were adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate. Genes with an adjusted P-value < 0.01 & Fold Change ≥ 2 found by DESeq2 were assigned as DEGs. GO enrichment analysis was implemented by the clusterProfiler packages based Wallenius non-central hyper-geometric distribution (Young et al., 2010). The clusterProfiler software were used to test the statistical enrichment of differential expression genes in KEGG pathways (Yu et al., 2012). Gene set enrichment analysis (GSEA) was performed with GSEA v3.0 (http://www.broadinstitute.org/gsea/).

To analyze the phylogenetic relationships of the accessions, we constructed an unrooted/rooted phylogenetic tree using the neighbor-joining method with the Kimura 2-parameter/p-distance model in MEGA-CC software (MEGAX) (Kumar et al., 2018), with 1000 bootstrap replicates. In total, 403,933 high-confidence SNPs were used to infer the population structure within accessions using ADMIXTURE (v1.22) (Alexander et al., 2009), K values (the putative number of populations) ranging from 1 to 10. We assessed the number of sub-populations using five-fold cross-validation. The Q matrix for each K value stacked assignment bar plots were generated using the R package Pophelper (http://royfrancis.github.io/pophelper). Principal component analysis (PCA) of the SNPs was performed using smartPCA program from the EIGENSOFT package using the default parameters (Price et al., 2006).

The SNPs markers were identified using the transcriptome data. Briefly, BAM files generated after the mapping process with HISAT2 were sorted using SAMtools (v0.1.19). Next, duplicate reads were removed with Picard tools (v2.25.7) (https://github.com/broadinstitute/picard). Finally, GATK (v3.8) was employed to call variants using the parameters “–indelSizeToEliminateInRefModel 50 –variant_index_type LINEAR –sample_ploidy 2 -nct 4 -U ALLOW_N_CIGAR_READS”. SNPs with a minor allele frequency (MAF) less than 0.01 or an integrity threshold below 0.5 were filtered out. Finally, 403,933 SNPs were retained for subsequent analysis. GWAS was performed using a Linear Mixed Model (LMM) in the GEMMA (V0.98.1) package (Zhou and Stephens, 2012). The matrix of pairwise genetic distances calculated by GEMMA and the PCA calculated by smartPCA were used as the variance-covariance matrix of random effects. The Manhattan and quantile quantile (QQ) plots of GWAS results were generated in R software (https://cran.r-project.org/).

To determine the relationship among candidate transcripts, eQTL analysis was performed by TRAS of 190 maize inbred lines using the LMM model. In our eQTL analysis, the 403,933 SNPs were defined as markers, and the expression of candidate transcripts were considered as phenotypes. To avoid potential false positives in multiple comparisons, the modified Bonferroni correction was used to determine the genome-wide significance thresholds of the TWAS, the more rigorous criterion of α=0.01/number of markers was used and the candidate genes were determined based on a threshold of P < 2.476e-08.

Total RNA was extracted from maize kernels using the FlaPure Plant Total RNA Extraction Kit following the manufacturer’s protocol. Subsequently, after DNase treatment, Reverse transcription was performed using the UnionScript First-strand cDNA Synthesis Mix. The resulting cDNA was then utilized as a template for quantitative PCR, which was performed with the GS AntiQ qPCR SYBR Green Fast Mix on the CFX96 Touch Real-Time PCR Detection System, including three biological replicates. The gene expression levels were normalized to the expression of theβ-actin and calculated for each sample using the 2^−ΔΔCT method. The information of the primers is shown in Supplementary Table S2 .

To investigate the genetic basis for folate characteristics in maize, we assembled a germplasm collection comprising 380 kernel samples derived from 190 maize inbred lines representing diverse geographical origins, including China, Thailand, Mexico, and other regions ( Figure 1A ). Transcriptome sequencing analysis were conducted on a total of 380 samples. A total of 2593.07 Gb of clean data reads with a Q30 score of 94.03% were generated, with average sequencing data per sample of 5.71 Gb. Population structure analysis based on transcriptome-derived SNP markers revealed distinct genetic characteristics among four regional subpopulations ( Figure 1B ). To investigate the genetic basis of folate accumulation, we quantified kernel folate content using two biological replicates per sample. Samples were stratified into high- and low-folate groups based on median folate levels, demonstrating significant inter-group divergence in folate content (P < 0.01, Figure 1C ). PCA of the folate-stratified groups showed consistent clustering patterns with phylogenetic relationships ( Figure 1D ), confirming that phenotypic variation in folate content occurred independently of genetic background ( Figure 1E ).

Using a threshold of Fold Change ≥ 2 and P value < 0.01, we identified DEGs between the high-folate and low-float groups, resulting in a total of 137 DEGs, of which 112 were upregulated and 25 were downregulated ( Figure 2A, B ). To elucidate the functional of DEGs, we performed GO enrichment analysis using a significance threshold of P < 0.05 across three categories: Molecular Function (MF), Biological Process (BP), and Cellular Component (CC) ( Figure 2C ). In the MF category, the top five enriched GO terms were: protein self-association, unfolded protein binding, acid phosphatase activity, L-3-cyanoalanine synthase activity, glucan endo-1,3-beta-glucanase activity, and phosphatase activity. These results suggest that the DEGs may participate in biochemical pathways associated with cyanide amino acid metabolism, phosphatase-related metabolic regulation, and protein homeostasis. In the BP category, the most significantly enriched terms included cyanide catabolic process, DNA catabolic process, protein folding, and response to reactive oxygen species (ROS). These findings indicate that the DEGs are likely involved in critical biological mechanisms such as DNA damage repair, stress responses to external stimuli, and metabolic regulation. In the CC category, the top five enriched GO terms were plasmodesma, lysosome, apoplast, extracellular space, and integral component of membrane. This implies that the DEGs may contribute to intercellular communication, substrate transport, and maintenance of cellular architecture. KEGG enrichment analysis revealed the DEGs were enriched in phenylpropanoid biosynthesis, minoacyl-tRNA biosynthesis, carotenoid biosynthesis, ceatin biosynthesis, and cyanoamino acid metabolism ( Figure 2D ).

To further elucidate the regulatory mechanisms of folate expression in maize endosperm, we conducted a GWAS using the LMM method implemented in GEMMA software and analyzed the candidate gene. We employed a significance threshold adjusted by the 1% Bonferroni correction (P-value ≤ 1.00E-05) to identify significant associations. By integrating the Manhattan plots for folate traits and LD decay rates of 10 chromosomes in 380 maizes, and based on the LD coefficient decreasing to half of its maximum at a distance of 100 kb, we selected target intervals at 2 kb upstream and downstream of the SNP, and finally identified 692 significantly associated loci ( Figure 3B ). These loci included 3,943 SNPs associated with folate traits ( Supplementary Table S3 ). Based on analysis of genes in LD regions, a total of 2,153 candidate genes were identified for potential folate characteristic associations.

We further identified the core genes associated with folate in maize by integrating GWAS with DEGs ( Figure 3A ). They were seven genes, including Zm00001eb252460, Zm00001eb049140, Zm00001eb242900, Zm00001eb277780, Zm00001eb402550, Zm00001eb317130, Zm00001eb253750, were found to be located on the SNPs identified in the GWAS ( Figure 3B ), and showed significant differences in expression levels between the low and high folate groups ( Figures 3C, E ). Therefore, these genes were identified as candidate genes associated with folate accumulation. We performed GSEA to further understand the functions of the hub gene, and the results revealed that these genes are significantly enriched in several biological pathways, including biosynthesis of amino acids, phenylpropanoid biosynthesis, oxidoreductase activity, and unfolded protein binding ( Figure 3D ).

To further investigate the causal genes associated with the seven core genes, we performed TWAS to identify potential genes linked to these core genes. In total, 13 candidate regulatory genes were detected with a significance threshold of P < 2.476e-08 ( Figure 4A ). It was worth noting that these genes identified by TWAS were significantly associated folate synthesis. Zm00001eb112740 and Zm00001eb169580 were annotated with dihydrofolate reductase activity and participation in the tetrahydrofolate (THF) biosynthetic process, indicating their central roles in folate metabolism by directly influencing the production of the active form of folate. The Zm00001eb107660 gene was linked to the biosynthesis of folate-containing compounds, 5-formyltetrahydrofolate cyclo-ligase activity, and tetrahydrofolate interconversion, suggesting its regulatory role in the synthesis and metabolic transformation of folate and its bioactive derivatives. Additionally, multiple genes, including Zm00001eb432940, Zm00001eb067370, Zm00001eb146170, Zm00001eb170020, Zm00001eb319990, and Zm00001eb404490, were implicated in THF synthesis and interconversion processes. Notably, Zm00001eb115340, Zm00001eb215470, Zm00001eb301040, and Zm00001eb285930 remain unannotated but represent potential candidates for folate biosynthesis, warranting further functional validation. The KEGG annotation analysis of the 13 identified genes is presented in Figure 4B ; Supplementary Table S4 .

To further analyze the correlation between these candidate genes and folate expression, we extracted RNA from 32 maize kernels each from both the high and low folate groups and performed qRT-PCR analysis to examine the expression levels of these candidate genes across the two groups. The results revealed that these seven core genes exhibited significantly higher expression levels in the high-folate group compared to the low-folate group ( Figure 5 ), consistent with the transcriptomic data analysis. Furthermore, we also conducted expression analysis on 13 causal genes selected from these core genes through qRT-PCR. The results showed that 11 out of 13 genes, including Zm00001eb107660, Zm00001eb112740, Zm00001eb115340, Zm00001eb146170, Zm00001eb169580, Zm00001eb215470, Zm00001eb285930, Zm00001eb301040, Zm00001eb319990, Zm00001eb404490, and Zm00001eb432940, were significantly higher in the high folate groups compared to the low folate group (P < 0.05). The remaining two genes exhibited differences between the two groups, but these differences were not statistically significant ( Supplementary Figure S1 ).