Time-resolved multiomics analysis of the genetic regulation of maize kernel moisture

Jianzhou Qu,Shutu Xu,*,Xiaonan Gou,Hao Zhang,Qian Cheng,Xiaoyue Wang,Chuang Ma,Jiquan Xue,*

a The Key Laboratory of Biology and Genetic Improvement of Maize in Arid Area of Northwest Region,Northwest A&F University,Yangling 712100,Shaanxi,China

b Maize Engineering & Technology Research Centre,College of Agronomy,Northwest A&F University,Yangling 712100,Shaanxi,China

c State Key Laboratory of Crop Stress Biology for Arid Areas,Center of Bioinformatics,College of Life Sciences,Northwest A&F University,Yangling 712100,Shaanxi,China

Keywords:Maize Kernel moisture Kernel dehydration rate GWAS Multiomics

ABSTRACT Maize kernel moisture content(KMC)at harvest greatly affects mechanical harvesting,transport and storage.KMC is correlated with kernel dehydration rate(KDR)before and after physiological maturity.KMC and KDR are complex traits governed by multiple quantitative trait loci(QTL).Their genetic architecture is incompletely understood.We used a multiomics integration approach with an association panel to identify genes influencing KMC and KDR.A genome-wide association study using time-series KMC data from 7 to 70 days after pollination and their transformed KDR data revealed respectively 98 and 279 loci significantly associated with KMC and KDR.Time-series transcriptome and proteome datasets were generated to construct KMC correlation networks,from which respectively 3111 and 759 module genes and proteins were identified as highly associated with KMC.Integrating multiomics analysis,several promising candidate genes for KMC and KDR,including Zm00001d047799 and Zm00001d035920,were identified.Further mutant experiments showed that Zm00001d047799,a gene encoding heat shock 70 kDa protein 5,reduced KMC in the late stage of kernel development.Our study provides resources for the identification of candidate genes influencing maize KMC and KDR,shedding light on the genetic architecture of dynamic changes in maize KMC.

1.Introduction

Maize(Zea mays L.)is a food and fodder crop worldwide.Kernel moisture content(KMC)affects maize production and breeding.High KMC at harvest increases the risk of ear sprouting,ear rot,and plant lodging,thus affecting the harvesting and safe storage of maize kernels,especially in short-to mid-season growing areas affected by environmental factors such as temperature,rainfall,and relative humidity[1–3].Reducing KMC at harvest has become a target of maize breeding and biotechnology-assisted improvement.

KMC at harvest is determined by two phases:the first phase is before physiological maturity,in which KMC is associated mainly with accumulation of dry matter via kernel filling and regulated mainly by genetic factors.The second phase is after physiological maturity,in which KMC is influenced by external environmental factors in addition to genetic factors[4,5].Kernel dehydration rate(KDR)also determines KMC at harvest,and maize varieties with fast KDR before and after physiological maturity generally have lower KMC at harvest[6,7].Both KMC and KDR are closely associated with kernel composition(including sugars,water-soluble polysaccharides,and some hydrophobic compounds)[8,9]and various agronomic traits(including kernel row number,ear length and diameter,pericarp thickness,and husk length and tightness)[10–15].Maize KMC and KDR are quantitative traits,making the identification of their genetic basis difficult.

In association studies aimed at identifying the genetic basis of KMC and KDR,multiple quantitative trait loci(QTL)for maize KMC and KDR have been detected in diverse populations by phenotype–genotype association analysis[1,3,16–28].However,no QTL have been fine-mapped or even cloned,owing to technical limitations and the inherent complexity of KMC and KDR.Further investigations are needed to elucidate the genetic basis of KMC and KDR and to perform gene-based breeding to improve maize KMC.

A genome-wide association study(GWAS),a high-resolution method,provides an opportunity to methodically investigate the genetic architecture of complex traits and identify beneficial alleles based on high diversity and rapid linkage disequilibrium(LD)decay[29,30].For maize KMC and KDR,several studies[28,31–33]have employed GWAS to identify favorable alleles.These studies showed the promise of linking GWAS loci to genes associated with KMC and KDR.With the development of next-generation sequencing technologies,individual omics data provide complementary support so that integrating multiomics data is expected to strengthen the signal for pinpointing genes associated with KMC and KDR.At a different omics level,multiomics data for individual genes are further amplified when multiple genes are considered together,particularly given the polygenicity of KMC and KDR.

The object of the present study was to identify gene sets influencing KMC and KDR by identifying candidate genes for KMC and KDR using large-scale GWAS and temporal transcriptomic and proteomic atlases.In promising candidate gene sets,we cloned and identified ZmHSP5,which encodes a heat shock 70 kDa protein 5.These results confirm and greatly expand our previous understanding of the genetic architecture of maize KMC and KDR and provide a resource for future research on the molecular mechanisms of reducing maize KMC at harvest.

2.Materials and methods

2.1.Plant materials and field design

An association mapping panel(named AM212)of 212 maize inbred lines comprised 197 inbred lines from the Shaan A(41)and Shaan B(156)groups[34]and 15 inbred lines collected domestically and abroad.Trials were performed at the maize breeding testing station of Northwest A&F University in Yangling,Shaanxi province(108°05′N,34°32′E)in 2018 and 2019.The field experiment followed a randomized complete block design with two replications.Each inbred line was planted in 4 rows with row lengths of 4.5 m and row spacing of 0.6 m.The planting density was 75,000 plants per hectare and standard irrigation and fertilization management practices were used throughout the developmental period.An ethyl methanesulfonate(EMS)mutant of Zm00001d047799(EMS4-0be358)was ordered from the Maize EMS induced Mutant Database(https://elabcaas.cn/memd/public/index.html#)[35],and was backcrossed once with B73 and then selfed once.The mutation site was genotyped by sequencing with primers listed in Table S1.

2.2.Phenotype data collection and analysis

The ears of all maize inbred lines were bagged before silk emergence.To ensure that all ears of AM212 were synchronized with respect to developmental stage to minimize any environmental influence on kernel dehydration,we performed self-pollination or sib pollination for three days and marked the date of pollination on the bag.KMC was measured at 10 successive stages(7,14,21,28,35,42,49,56,63 and 70 DAP).Three well-pollinated ears pollinated on the same day were collected in each plot and threshed by hand immediately after collection.One hundred kernels were excised from each ear with a scalpel and placed in an aluminum box(80 mm diameter×50 mm height).After sampling,kernel fresh weight(W1)was measured with a 0.001 g digital scale and each box was labeled with the date,sample name,and number and then placed in an oven.The samples were heated at 105 °C for 30 min and then dried at 70 °C to constant weight(W2).KMC was calculated as.

Two methods were used to calculate KDR:KDR1,

where n is the nth sampling date,and KDR2,the area under the dry down curve(AUDDC):

where n is the number of assessments,γ is the mean KMC,i is the ith sample,and t is the pollination date(7,14,21,28,35,42,49,56,63 and 70 DAP)[36].

A mixed linear model was built by fitting intercept as a fixed effect and genotype and year as random effects in SAS 8.1(SAS Inc.,Cary,NC,USA).The best linear unbiased prediction(BLUP)value for each maize inbred line was calculated.Broad-sense heritability(H2)was calculated as follows:

where δ2Gis genotypic variance;δ2Eis residual variance;and n is the number of years[37].All data from all developmental time points of each year and BLUP were used for subsequent analyses,including 30 KMC,135 KDR1,and 135 KDR2 traits.

2.3.DNA sample preparation,genotyping,and filtering

Fresh young leaves were collected 3 weeks after sowing and quickly frozen in liquid nitrogen.Total DNA was manually extracted by the cetyltrimethylammonium bromide method[38].Total DNA was quantified with a NanoDrop ND-2000(Thermo Fisher Scientific,Waltham,MA,USA)and 1% agarose gel electrophoresis.DNA samples that passed the quality check and then performed genotype detection using the Affymetrix maize 6H90K Chip(Beijing Compass Biotechnology Co.,Ltd.,Beijing,China).A total of 73,006 markers were generated,and SNPs were filtered using a missing rate≥20% and minor allele frequency(MAF)≤5% using PLINK[39].Missing genotypes were imputed with Beagle 5.1[40]with default parameters.

2.4.Genome-wide association analysis

GWAS was performed using a mixed linear model implemented in GEMMA 0.98.4[41].PCA was performed using BLUP values with the FactoMineR and factoextra packages[42].KMC,KDR1,and KDR2 data from all years,BLUP data,and the first 10 PCs for all traits were subjected to GWAS.The number of PCs of the associated population was estimated with GCTA[43].A neighbor-joining tree was constructed using MEGA-X 10.0.5[44]and visualized with iTOL(https://itol.embl.de).The Bonferroni-corrected threshold(0.05/the number of effective SNPs)was shown to be too stringent[45].Accordingly,a P-value of 1.49×10-5(P=1/N,N=67,076)was adopted as a threshold for declaring a significant association by the adjusted Bonferroni method[46].The mean LD decay distance was calculated using PLINK 1.90b6.25 and a value of r2＜0.2 was considered to indicate linkage.The mean r2for all chromosomes was estimated at～200 kb,when the value of the cutoff for r2was set to 0.1.For each significant SNP detected,a 200-kb region flanking the SNP was searched for candidate genes.The physical location of the SNPs was identified based on the maize genomic sequence version 3.0(ZmB73_RefGen_v3;https://www.maizesequence.org),and all v3 candidate gene IDs were converted to v4 gene IDs by the MaizeGDB database(https://www.maizegdb.org).

2.5.RNA-seq and data analysis

Eight maize inbred lines were selected from fast-dehydration types(KB182,KA225,KA105,and KB035,designated FD1–FD4)and slow-dehydration types(PH6WC,KB020,KA327,and 91227,designated SD1–SD4).The kernels of these inbred lines were sampled by manual dissection at 21,28,35,42,49,and 56 DAP in 2018,frozen immediately in liquid nitrogen,and stored at-80°C.Three biological replicates were set up for each developmental time point.Total RNA was extracted with TRIzol reagent(Invitrogen,Carlsbad,CA,USA)following the manufacturer’s instructions.RNA-seq libraries were constructed according to the manufacturer’s protocol of the Illumina NEBNext Ultra RNA Library Prep Kit(Illumina,Inc.,San Diego,CA,USA).RNA sequencing with three biological replications using an Illumina NovaSeq-PE150 Platform(Illumina,Inc.),which was completed by Novogene(Novogene Co.,Ltd.,Tianjin,China).Sequencing reads were quality trimmed with Trimmomatic 0.38[47].After low-quality read removal,the remaining reads were aligned to the maize B73 reference genome(RefGen_v4; ftp://ftp.ensemblgenomes.org/pub/plants/release-41/fasta/zea_mays/dna)sequence assembly with HISAT2 2.0.5[48].StringTie 1.3.5(https://ccb.jhu.edu/software/stringtie)was used to normalize and estimate gene expression levels in fragments per kilobase of transcript per million mapped reads(FPKM).The FPKM values of triplicate samples were averaged for each gene.Differentially expressed genes(DEGs)were identified based on false discovery rate(FDR)-adjusted P-value≤0.05 and a fold change≥2 or≤0.5 using edgeR[49].Enrichment analysis of Gene Ontology(GO)terms was performed with AgriGO 2.0(https://bioinfo.cau.edu.cn/agriGO/analysis.php)[50].Transcription factors(TFs)were identified by reference to GRASSIUS(https://grassius.org)[51]and PlantTFDB(https://planttfdb.cbi.pku.edu.cn)[52].

2.6.Identification and quantitation of proteins

Maize kernels of maize inbred lines KA105(FD3)and 91,227(SD4)were sampled at 28,35,42,49,and 56 DAP.Protein extraction and peptide preparation were performed by Novogene(Novogene Co.,Ltd.).Three biological replicates of the samples at each time point were mixed as a pooled sample,and this pooled sample was used to construct the library for the data-independent acquisition(DIA)protein identification and the Mass Spectrum(MS)analysis for each sample as described in a previous study[53].

The data-dependent acquisition(DDA)MS raw files were analyzed with Proteome Discoverer 2.2,and peak lists were searched against the SwissProt database(https://www.ebi.ac.uk/UniProt).MS1-based label-free quantification(LFQ)was performed using the maxLFQ algorithm[54].MS2-based label-free quantification was performed by analyzing DIA raw data using Biognosys Spectronaut 9.0.DIA MS/MS quantification is the sum of the integrated fragment ion peak areas.Subsequently,data were further analyzed as described in a previous study with minor modifications[55].At least one unique peptide with an FDR≤1% was required for protein identification and quantification.In the comparison group,a protein was assigned as a differently expressed protein(DEP)using the cutoff criteria of an FDR-adjusted P-value≤0.05 and a fold change≥1.5 or≤0.67.To confirm quantitative protein,20 proteins were randomly chosen and quantified by PRM analysis at Novogene BioLabs(Novogene Co.,Ltd.).

2.7.WGCNA

WGCNA was performed by the WGCNA R package(version 1.63; https://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA).The DEGs and DEPs identified across developmental time points were used for module construction.All gene and protein abundances were normalized.Pearson’s correlation coefficients(PCCs)for each gene–gene and protein–protein comparison were calculated,and an adjacency matrix of the connection strengths was constructed.The best power β was optimized to adjust the scale-free property of the coexpression network and the sparsity of connections between genes or proteins.Highly similar clusters were merged in the network using the mergeCloseModules function using a cutHeight value of 0.15.Module–trait associations were estimated using the correlation between the module eigengene and the phenotype(|PCC|≥0.7,P-value≤0.01),which allows easy identification of the expression set(module)highly associated with the phenotype.Hierarchical clustering was performed with the R package pheatmap(https://cran.r-project.org/web/packages/pheatmap/index.html)using PCC as the distance measure.Venn diagram samples were visualized with the R package UpSetR(https://cran.rproject.org/web/packages/UpSetR/index.html).The correlation relationship between samples was calculated with the R package Hmisc (https://cran.r-project.org/web/packages/Hmisc/index.html).

3.Results

3.1.Phenotypic variation in KMC and KDR

We observed abundant variation in KMC and KDR in the AM212 panel.The inbred lines showed substantial variation in KMC,which followed a near-normal distribution,with a range of values of 4.61%–31.66% at 10 developmental time points in the two years(Fig.S1;Table S2).The KMC values were strongly correlated between the two years(Fig.S1A).The KMC data of 10 developmental time points were further transformed to KDR1 and KDR2 values using two methods(Fig.1A).Analysis of variance suggested that the variation of maize KMC,KDR1 and KDR2 was influenced mainly by developmental stage,followed by genotype,environmental factors,and additional interaction effects(Fig.S1B).KMC,KDR1 and KDR2 showed high broad-sense heritability(H2),ranging from 70.67% to 85.53% for KMC,14.77% to 80.82% for KDR1,and 79.03% to 89.79% for KDR2(Fig.S1C).

We estimated the BLUP value of each inbred line across years at each time point for KMC,KDR1,and KDR2.The BLUPs showed normal distribution and were consistent with the variation in KMC,KDR1,and KDR2 at two years(Figs.1B,S2).The range of KMC variation in AM212 continuously decreased from the early to late stages,whereas the KDR at the early stage of kernel development was faster than that at the late stage(Fig.1C).Cluster analysis revealed that AM212 KMC variation could be broadly divided into three types:inbred lines with slow dehydration(type I),progressive dehydration(type II),and fast dehydration(type III)(Fig.S3).

Additional phenotypic data were generated from the time series phenotypic data(KMC,KDR1 and KDR2)for principal component analysis(PCA)-based GWAS analysis,inspired by previous research[56].For KMC,the first two principal components(PCs)explained respectively 63.32% and 12.32% of the total variance of development points,and the KMC at 42 DAP showed high positive loadings on PC1(78.08%)(Fig.1D),suggesting that 42 DAP may be a key time point for KMC transformation in late kernel development.In contrast,the first 10 PCs of KDR1(99.81%)and KDR2(99.92%)explained nearly 100% of the KDR variance among developmental stages(Figs.1D,S4),and all PCs displayed roughly normal distributions(Fig.1E).Accordingly,these first 10 PCs of KMC,KDR1,and KDR2 were used as quantitative traits for GWAS.

Fig.1.Phenotypic diversity of 212 inbred lines.(A)Collection procedure of phenotype data.(B)Distribution of BLUP values for KMC and KDR.(C)Changes in BLUP values of KMC and KDR during kernel development.The change in KMC in individual inbred lines is depicted in gray lines,and colored lines indicate the overall pattern after fitting.(D)Loading plot of PC1 and PC2 based on BLUP values of KMC and KDR.Purple,green and yellow indicate clusters of developmental stages.The proportions of variance for PC1 and PC2 are shown in parentheses.(E)Distribution of PC scores of KMC and KDR.KDR was calculated by two methods and named KDR1 and KDR2;see details in the Methods section.

3.2.Genetic architecture of KMC and KDR

To characterize the genetic background of the association panel,we generated a final set of 67,076 high-quality single nucleotide polymorphisms(SNPs)from the 6H90K Chip with a missing rate＜0.20 and MAF＞0.05.PCA of these 67,076 SNPs showed that maize inbred lines from the Shaan A group and the Shaan B group were mostly well separated,and the 15 introduced inbred lines clustered with the Shaan B group(Fig.S5A).This observation was also reflected by the phylogenetic relationships of the 212 inbred lines,which were further divided into five main clades(Fig.S5B).The LD decay distance was approximately 200 kb(r2=0.20)(Fig.S6),in agreement with a previous study using other maize populations[57].

To investigate the genetic basis of KMC and KDR,we performed a GWAS using genome-wide efficient mixed-model association(GEMMA)software for KMC and KDR,including all years and BLUP and PCA data(Fig.2A–C).In total,98 significantly associated loci(-log10P≥4.83)for KMC were detected with phenotypic variation explained(PVE)ranging from 0.22%to 22.92%(Fig.2D–E),of which 20 loci were consistently detected at multiple time points,and approximately 44% of KMC-associated loci(43/98)were located in previously reported QTL intervals(Table S3).We detected 222 loci for KDR1 with 0.09%to 22.90%PVE and 107 loci for KDR2 with 0.71% to 23.51% PVE(Fig.2D–E).Among these associated loci,97 loci were repeatedly detected in multiple GWAS results,and 141 loci were located in previously reported QTL intervals(Table S3).Specifically,33 and 214 loci were uniquely detected for KMC and KDR(Fig.2D;Table S3),indicating that the variation in KMC and KDR was also associated with differences in their genetic basis.Of 85 loci,70 were detected using PCA-transformed phenotypic data instead of the BLUP values of KMC and KDR(Table S3).

Fig.2.GWAS for KMC and KDR.(A–C)The colors of the dots correspond to a single trait,including all yearly and BLUP and PCA data.SNPs above red lines exceed a significance threshold(-log10P≥4.83).Only SNPs with-log10P≥2 are plotted.Manhattan plot of GWAS associations for KMC(A),KDR1(B)and KDR2(C).(D)Colocalization of SNPs detected among KMC,KDR1 and KDR2.(E)Distribution of effect values of significant SNPs for KMC,KDR1,KDR2 and overall KDR.(F)Numbers of potential candidate genes for KMC and KDR traits.

We further calculated the additive effects of significant loci for KMC and KDR and found that these significant loci had 3.21%–42.70% additive effects for variance of KMC and 2.79%–42.82%additive effects for variance of KDR.The mean additive variance of significant loci for KDR was higher than that for KMC,far lower than the broad-sense heritability(Fig.S7).These findings suggest that KMC and KDR are controlled by a few large-effect and many small-effect additive SNPs.We estimated a candidate region from the flanking 200-kb regions of each significant SNP for KMC and KDR.These regions contained 197 unique potential candidate genes for KMC and 1363 unique potential candidate genes for KDR,with an overlap of 446 genes(Fig.2F;Table S4).

3.3.Identification of module genes and proteins highly associated with KMC

To systematically investigate the relationship between gene/protein activity and dynamic KMC,we performed highthroughput RNA-Seq and whole-proteome profiling for the kernels of eight maize inbred lines with diverse dehydration types(Figs.3A,S8).At the transcription level,26,393 genes were expressed with a mean FPKM≥1 of three biological replicates at a minimum of one time point.After removing genes without significant expression changes(FDR-adjusted P-value＞0.05 and/or|log2(fold change)|＜1),we retained 17,435 genes with differential expression between time points in the same inbred line for the downstream weighted gene coexpression network analysis(WGCNA)with the R package WGCNA[58](Fig.3A;Table S5).

The KMC-associated modules were identified based on a Pearson correlation coefficient(PCC)≥0.7 and a P-value≤0.001.To further select promising genes for KMC,we considered only genes that were uniformly detected in at least three FD or SD inbred lines,and a PCC＞0.7 with KMC was used for subsequent analysis(Fig.3A).As a result,we identified 1406 and 608 genes that were unique to FD and SD,respectively,and 1097 genes overlapping between FD and SD,including 92 TFs(Fig.3B;Table S6).GO enrichment analysis showed that KMC-associated module genes detected specifically for FD were enriched mainly in starch biosynthetic process,stomatal movement,metal ion transport,and response to light stimulus(Fig.3C;Table S7).In contrast,KMC-associated module genes detected specifically for SD were significantly enriched in brassinosteroid homeostasis,lipid metabolic process,and coumarin biosynthetic process(Fig.3C;Table S7).We identified genes common to FD and SD that function in kernel development,including seed dormancy process and organ senescence;steroid and alcohol biosynthetic processes;chlorophyll catabolic process and lipid storage;nitrate and water transport;and multiple stress responses,including freezing,heat acclimation,salt stress,desiccation,hydrogen peroxide,abscisic acid(ABA),and fructose responses(Fig.3C;Table S7).

At the translation level,we applied both DIA and DDA workflows for whole-proteome profiling that quantified 6482 gene productions with a 1% peptide FDR in at least two of three replicates per time point(Fig.3A).A total of 3700 differentially expressed proteins(DEPs,P-value≤0.05 and fold change≥1.5 or≤0.67)were identified by pairwise comparisons,similar to those used in transcriptome comparison(Table S8).We performed PRM analysis on 20 proteins,revealing expression transitions similar to the DIA results(Fig.S9).

Fig.3.KMC-associated module genes and proteins and their functional diversity.(A)Workflow for identified module genes and proteins highly correlated with KMC.First,a coexpression network was constructed for each maize inbred line using temporal transcriptome and proteome data.Second,modules highly correlated with KMC were identified(PCC≥0.70 and P-value≤0.01).Third,module genes and proteins highly correlated with KMC were selected(|PCC|≥0.70 and P-value≤0.01).Finally,comparison of KMC-associated module genes and proteins in FD and SD.(B)Comparison of KMC-associated module genes in FD and SD.Red numbers in parentheses indicate TFs.(C)GO enrichment for KMC-associated module genes.(D)Comparison of KMC-associated module proteins in FD and SD.Red numbers in parentheses indicate TFs.(E)GO enrichment for KMC-associated module proteins.

By comparing the proteins associated with KMC between FD and SD,we found 559(including 10 TFs)and 155 unique proteins(including 1 TF)for FD and SD,respectively,and 45 proteins(including 1 TF)overlapping between FD and SD(Fig.3D;Table S9).GO enrichment analysis showed that unique proteins associated with FD play roles in kernel development,including proline,short-chain fatty acid and starch biosynthetic processes;water transport,callose deposition in the cell wall,cell wall pectin metabolic process,NADPH regeneration,the tricarboxylic acid cycle and gluconeogenesis;and response to heat,desiccation,hydrogen peroxide and cytokinin(Fig.3D;Table S10).In contrast,unique proteins for SD were linked to lipid localization,killing of cells of other organisms and response to hypoxia and oxidative stress(Fig.3D;Table S10).Overlapping proteins were associated with cellular amine and sulfur compound metabolic processes and hormone and pigment biosynthetic processes(Fig.3D;Table S10).

3.4.Identification and validation of candidate genes associated with KMC and KDR

Integrating GWAS analysis,network-based temporal transcriptome data analysis and network-based temporal proteome data analysis, we identified four genes (Zm00001d002258,Zm00001d002260,Zm00001d027290,and Zm00001d035920)that were consistently linked to KMC and KDR at the multiomics level(Fig.4A).Zm00001d035920,encoding an RNA binding protein,is a binding partner of ACD11.1,hereafter named ZmBPACD11.Within this gene,two adjacent SNPs(Affx-291435855 and Affx-291445022)were significantly associated with KMC and KDR(Fig.4B–C).ZmBPACD11 showed lower gene expression levels but higher protein abundance in FD than in SD inbred lines(Fig.4D).Protein–protein interaction analysis showed that ZmBPACD11(B4FEG8)has potential interactions with six proteins,of which A0A1D6DZ56 encodes an ABC transporter and has been verified to be involved in diverse processes,such as pathogen response,surface lipid deposition,phytate accumulation in seeds,and transport of the phytohormones auxin and ABA(Fig.4E)[59].This result suggests that ZmBPACD11 is likely to be a key candidate gene for KMC and KDR.

Fig.4.ZmBPACD11 is a potential candidate gene for KMC and KDR.(A)UpSet plot showing the number of genes associated to KMC and KDR that were detected by multiomics.(B)The upper panel indicates GWAS of the BLUP value for KMC at 56 DAP and KDR2 from 49 DAP to 56 DAP.The lower panel indicates the P values of SNPs surrounding a 1-Mb region upstream and downstream of the lead SNP on chromosome 6,gene model of ZmBPACD11 and LD pattern of the ZmBPACD11 region.The red stars mark the positions of the lead SNP.(C)The influence of the lead SNPs on KMC and KDR.The significance of the difference between two genotypes was evaluated by two-tailed Student’s t-test.(D)Expression patterns of ZmBPACD11 at the transcriptional and translational levels.(E)Protein interaction network of ZmBPACD11.

Another GWAS locus on chromosome 9 was associated with KDR at multiple kernel developmental stages.In the candidate region,there was one candidate gene Zm00001d047799,transcriptional and translational changes of which were significantly associated with changes in KMC(Fig.5A–C;Tables S6,S9).Zm00001d047799 encodes heat shock 70 kDa protein 5,which has been reported to function in protein folding,calcium homeostasis,and serves as an essential regulator of the endoplasmic reticulum(ER)stress response and autophagy[60,61]and is hereafter named ZmHSP5.This gene showed a higher gene expression level and lower protein abundance at the late stage of kernel development in FD than in SD(Fig.5D).To further determine the potential function of ZmHSP5,we identified an EMS mutant(EMS4-0be358)carrying a C-to-T substitution in the second exon that results in CAG becoming the stop codon TAG,thus introducing a premature stop codon(Fig.5E).KMC was significantly lower in homozygous HSP5 plants than in their homozygous wild-type siblings,and KDR was also faster in homozygous HSP5 plants than in their homozygous wild-type siblings at the late stage of kernel development(Fig.5F).WGCNA showed that 120 genes were consistently detected in the coexpression networks of at least two maize inbred lines and synchronized with ZmHSP5 expression,and they were enriched mainly in regulating kernel development and responding to multiple stresses,including seed dormancy process and response to ABA,hydrogen peroxide,desiccation,salt stress,cold,heat acclimation,and high light intensity(Fig.S10;Table S11).These results indicate that the ZmHSP5 function affects KMC and,in turn,KDR.

4.Discussion

In this study,multiple highly promising candidate genes associated with KMC and KDR were identified,and these genes may drive multiple metabolic processes in maize kernel development.This large collection of genes provides a rich resource for investigating the genetic basis of KMC and KDR throughout maize seed development,cloning KMC-and KDR-associated genes,and designing precision breeding strategies.These data also provide many missing details for studying the spatiotemporal regulation of development processes in maize kernels.

Fig.5.ZmHSP5 affected KDR.(A)Manhattan plot of the GWAS for KDR2 from 14 DAP to 56 DAP in 2019.The lead SNP is indicated by a red dot.(B)The upper panel shows the P values of SNPs surrounding a 1-Mb region upstream and downstream of the lead SNP.The middle panel shows the gene model of ZmHSP5(v3).The lower panel shows LD pattern of the ZmHSP5 region.Red stars mark the positions of the lead SNP.(C)The effect of the genotype on KDR2 in the association-mapping population.(D)Expression patterns of ZmHSP5 at the transcriptional and translational levels.(E)Genotype of the wild type(WT)and HSP5 and the position of mutation on ZmHSP5(v4).The blue shading indicates the mutated nucleotide.In(C)and(E),the significance of the difference between two genotypes was evaluated by two-tailed Student’s t-test.(F)Quantification of KMC and KDR between the WT and HSP5 from 35 DAP to 49 DAP.

KMC and KDR are both complex quantitative traits regulated by multiple genetic factors and perturbed by a variety of environmental factors,including humidity,temperature,radiation,wind speed,and rainfall[4].Thus,obtaining an accurate phenotype is a prerequisite for KMC-and KDR-associated research.To eliminate the influence of environmental factors,the following measures were taken in this study:(1)inbred lines with similar growth periods were selected,(2)pollination dates were adjusted according to the flowering time of materials and centralized pollination,and(3)uniform samples were collected with the same pollination date.In addition,the selection of the measurement method is crucial for obtaining accurate KMC values,and many methods have been reported for measuring KMC,including the SK300 moisturedetermination meter[62],the digital timber-moisture meter and the electrophysics moisture meter model MT808[28,31,36].Although these methods are effective,it is necessary to consider the effective range of measurement,the position at which the probe penetrates the kernel,and whether precise calibration is required.We chose an improved oven-drying method with relatively high accuracy to measure KMC,although it is timeconsuming and labor-intensive[19,32,63].

In addition to the determined phenotype data,PCA is also an effective evaluation means for complex quantitative traits.PCA is valuable for extracting underlying factors for multiple traits by reducing dimensionality,and PC scores as dependent variables are proposed as a strategy for performing efficient GWAS[56,64–66].PC scores produced by PCA can transform skewed original variables into approximate normal distributions,resulting in robust,reliable GWAS results[67,68];GWAS using PC scores may detect genomic regions that could be overlooked by use of individual traits,given that PC scores represent integrated variables.In this study,we performed a GWAS and identified 98 and 279 loci associated with KMC and KDR,respectively(Table S3).Of these loci,a GWAS of 10 PCA datasets detected 15 loci consistent with the results of the GWAS for the yearly and BLUP datasets and identified 70 loci controlled by genetic factors that were detected specifically from PCA-transformed phenotypic data.

Most loci had small effects on KMC and KDR variation,suggesting that the sum of many SNPs with small effects is responsible for the main contribution to the genetic basis of the variations in KMC and KDR.In comparison with previous studies[16,19,27,62,69–72],approximately half of the loci were located in reported QTL intervals(Table S3).This finding could be the result of differences in genetic backgrounds,population size,captured recombinant events,environmental effects,and GWAS analytical methods.These loci identified during dynamic kernel development constitute a resource that can be used for accelerating mining of superior alleles for KMC and KDR.

Seed dehydration from the phase of reserve accumulation to maturation drying is associated with distinct gene/protein expression and metabolic switches[73].One critical question is which specific genes function in regulatory metabolic pathways and the gene regulatory network associated with maize KMC and KDR.Previous genetic studies have demonstrated the functional significance of certain genes in seed maturation and dehydration.Overexpressing the LEA Rab28 gene increased the water-stress tolerance of transgenic maize plants[74,75],and abi3 and vp1 mutations increased the desiccation tolerance of seeds[76,77].ATEM6 mutation and overexpression of AtMYB118 caused premature dehydration of seeds[78,79],and GAR2 affected maize KDR[28].However,such well-studied genes represent a small proportion of the genes associated with KMC and KDR that we have detected.We integrated GWAS,transcriptomics,and proteomics data and identified additional highly promising genes,such as ZmBPACD11 and ZmHSP5.This was especially true for ZmHSP5,which was confirmed by mutation.

The water loss of maize kernels is associated not only with the accumulation of storage material processes before physiological maturity.As water concentration declines,damaging processes can occur,and seeds must have,or acquire,mechanisms to protect against them.Seeds entering the late maturation phase have already gained their germination capacity and desiccation tolerance.Some biological processes have been found[73,80,81]to function in seed development and influence seed dehydration,including the accumulation of disaccharides,oligosaccharides,and the plant hormone ABA;synthesis of storage proteins,late embryogenesis-abundant proteins and heat-shock proteins;activation of antioxidative defenses;changes in the physical structure of the cell;and a gradual and steady increase in density.In this study,gene sets associated with KMC and KDR were enriched in a variety of unreported biological processes,including the biosynthesis process of fatty acids,glucan,starch and jasmonic acid;cell cycle;and glycolytic process,and they were also linked to various kernel desiccation tolerance processes,including response to water deprivation,heat,hydrogen peroxide,reactive oxygen species,fatty acid beta-oxidation and oxidation–reduction processes(Tables S6,S9).The content of amylopectin with a degree of polymerization of 6–12 was lower in inbred lines with fast KDR[82].These results provide a new perspective for understanding the biological processes that affect maize KMC and KDR.We provide detailed transcriptional and translational information on the spatiotemporal patterns of all these genes associated with KMC and KDR.Although the underlying molecular mechanism of these gene sets associated with KMC and KDR should be further investigated,multiple GWAS data combining spatiotemporal transcriptome and proteomic data provide genetic resources for mining key genes for KMC and KDR.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

CRediT authorship contribution statement

Jianzhou Qu:Investigation,Data curation,Methodology,Software,Visualization,Writing–Original draft.Shutu Xu:Conceptualization,Supervision,Validation,Writing–review & editing.Xiaonan Gou:Investigation,Validation.Hao Zhang:Investigation,Validation.Qian Cheng:Software,Validation.Xiaoyue Wang:Investigation,Validation.Chuang Ma:Conceptualization,Supervision,Visualization,Validation,Writing–review & editing.Jiquan Xue:Conceptualization,Supervision,Validation,Writing–review& editing.

Data availability

All phenotype data of inbred maize lines in this study are included in Table S1.All RNA-seq data generated in this study have been submitted to the NCBI Gene Expression Omnibus(GEO;https://www.ncbi.nlm.nih.gov/geo)under accession number GSE158816.The MS raw data generated in this study have been submitted to ProteomeXchange database(www.proteomexchange.org)via the iProX partner repository under accession number PXD021680.All other reasonable requests for data and research materials should be addressed to the corresponding authors.

Acknowledgments

This research was supported by Natural Science Foundation of Shaanxi Province(S2021-JC-WT-006),the National Key Research and Development Program of China(2018YFD0100200),the China Postdoctoral Science Foundation(2018M633588)and the China Agriculture Research System(CARS-02-77).We thank all members in the Key Laboratory of Biology and Genetic Improvement of Maize in Arid Area of Northwest Region for their helps and supports in field material planting and phenotype collection.

Appendix A.Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2022.04.017.