ILa kernels at 10 DAP had lower IAA levels and higher SA levels than ILb kernels
Two maize lines, ILa (JN14–7-22) and ILb (JN14–7-13), were developed from a single F2 progeny ear of a self-crossed hybrid after another three rounds of self-crossing. The mature dry kernel size of ILa selfing progenies was smaller than that of ILb kernels. The kernel length and width of ILa were smaller than those of ILb kernels. The one hundred-kernel-weight (HKW) of ILa kernels was significantly lower than that of ILb; the HKW of ILb was 21.56% more than that of ILa (Fig. 1). While the mature plant heights of ILa and ILb were similar (Additional file 2: Figure S1).
To investigate the possible role of phytohormones in determining kernel size during early kernel development in maize, the levels of the most important growth and defense hormones (IAA, SA, and JA) were measured in ILa and ILb kernels at 5, 8, and 10 DAP. At 5 DAP, no IAA was detected in the kernels of either inbred. At 8 DAP, IAA was detected in low abundance in kernels of both inbreds; the IAA level was slightly higher in ILa kernels than in ILb kernels. At 10 DAP, the IAA level was significantly lower in ILa kernels than in ILb kernels. SA was detected at relatively similar high abundance in ILa and ILb kernels at 5 and 8 DAP; at 10 DAP, SA was significantly higher in ILa kernels than in ILb kernels due to a sharp reduction of SA in ILb kernels and a relatively small reduction of SA in ILa kernels. JA was significantly lower in ILa kernels than in ILb kernels at 5 and 8 DAP, it was stable in developing ILa kernels at all three timepoints (Fig. 2). These results indicate that the significantly lower IAA level and higher SA level at 10 DAP in ILa kernels may correlate with the smaller kernel size and HKW of ILa kernels, compared to that in ILb kernels.
ILa and ILb kernel transcriptomes significantly differed at 5 and 8 DAP
We investigated the gene expression networks during early kernel development in ILa and ILb to evaluate potential molecular mechanisms regulating phytohormone levels and mature kernel sizes. Developing kernels were collected at 5, 8, and 10 DAP from field-grown ILa and ILb plants for comparative transcriptome profiling by performing principal component analysis (PCA). Striking differences were observed in ILa and ILb kernel transcriptomes at 5 and 8 DAP. The ILa kernel transcriptome at 5 DAP (a5) clustered closely to the ILb kernel transcriptome at 8 DAP (b8), suggesting that the temporal developmental progression was faster in ILa kernels than in ILb kernels. The ILa kernel transcriptome at 8 DAP (a8) and the ILb kernel transcriptome at 5 DAP (b5) were strikingly different from the other kernel transcriptomes. By contrast, the ILa kernel transcriptome at 10 DAP (a10) clustered closely to the ILb kernels at 10 DAP (b10), indicating that the difference between ILa and ILb kernels at 10 DAP was reduced (Additional file 3: Figure S2).
The parameters used for screening differentially expressed genes (DEGs) between ILa and ILb were fold change (FC) of the expression level in ILa [FC ≥2 or FC ≤0.5 under P-value ≤0.05, false discovery rate (FDR) ≤0.05)] compared to the expression level in ILb at the same developmental stage based on the sequenced fragments per kilobase of transcript per million mapped reads (FPKM). The transcriptional differences between ILa and ILb peaked at 8 DAP, as 1261 DEGs were found in the a8/b8 transcriptome pair, and 957 (76%) of these DEGs were upregulated in ILa kernels, compared to the expression levels in ILb kernels at 8 DAP. And 847 DEGs and 530 DEGs were obtained from the a5/b5 transcriptome pair and the a10/b10 transcriptome pair, respectively. Only 100 DEGs from these transcriptome pairs overlapped among all developmental stages; a5/b5 shared 292 and 159 DEGs with a8/b8 and a10/b10 transcriptome pairs, respectively. Only 165 DEGs overlapped between a8/b8 and a10/b10 transcriptome pairs (Fig. 3a, b). Of the overlapped 100 DEGs, the FPKM value of 47 DEGs were higher in ILa kernels at almost all time-points, including 16 DEGs whose FPKM value decreased to below 1 in ILb kernels after 8-DAP, and 35 DEGs were exclusively expressed in ILa kernels and 6 DEGs were exclusively expressed in ILb kernels at all time-points (Additional file 4: Figure S3). This indicates that almost all the overlapped 100 DEGs displayed similar direction of gene expression changes among different contrasts, except a few DEGs that showed opposite gene expression changes at 8- DAP or at 10-DAP. Gene ontology (GO) analysis indicated that most of these DEGs were enriched in biological processes related to responses to abiotic stress (temperature, heat, light, water, osmotic stress, chemical stimulus, salt, auxin, ABA, and ROS) or biotic stress (defense responses to fungi and bacteria). The enrichment of defense-related functional categories was consistent with the high levels of SA and JA in maize kernels at 5 and 8 DAP (Fig. 2). The top-enriched molecular functions were catalytic activity, hydrolase activity, and oxidoreductase activity (Fig. 3c). The a5/b5 and a8/b8 transcriptome pairs contained 54 and 68 DEGs, respectively, that encode components of ‘response to hormone stimulus’, indicating significant differences in hormone signaling in ILa and ILb kernels during this period.
The top-enriched KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways of these DEGs were secondary metabolite biosynthesis pathways, including flavonoid, phenylpropanoid, flavone and flavonol, stilbenoid, diarylheptanoid and gingerol, and benzoxazinoid biosynthesis pathways. Other highly enriched pathways were metabolic pathways, including ascorbate and aldarate metabolism, starch and sucrose metabolism, and phenylalanine metabolism. Large number of DEG were enriched in ‘plant-pathogen interaction’ pathways at all timepoints. Abundant DEGs in a5/b5 and a8/b8 transcriptome pairs were enriched in ‘starch and sucrose metabolism’, ‘plant hormone signal transduction’; by contrast, these two pathways were remarkably reduced in DEG numbers in a10/b10 transcriptome pairs (Fig. 3d). Of the 61 DEGs that enriched in ‘plant hormone signal transduction’ from a8/b8, most of them expressed higher in ILa than that of in ILb at almost all the three time-points, including 7 of them specifically expressed in ILa kernels; except 17 of them expressed higher in ILb than that in ILa at almost all the time-points (Additional file 5: Figure S4). These results indicate that there are remarkable differences in hormone- and defense-related GO/KEGG functional categories in ILa and ILb kernels at 5 and 8 DAP, suggesting that phytohormone signaling and defense responses may function as molecular regulators mediating early maize kernel development.
Differences in phytohormone signaling and defense response in ILa and ILb kernels peaked at 8 DAP
Consistent with the above results, ILa and ILb kernels at 8 DAP had striking differences in the expression of certain transcription factor (TF) encoding genes, auxin-, kernel-development- and defense- related genes (Fig. 4). We identified 102 TF-encoding DEGs that displayed significantly higher expression levels in ILa transcriptomes than in the respective ILb transcriptomes, especially at 8 DAP. The TF-encoding DEGs were from 18 families, including MYB, MADS, bZIP, NAC, AP2/EREBP, WRKY, and GATA. The Yuc1 expression level was significantly higher in a8 than in b8 kernels, whereas Yucca2 expression was higher in b8 than in a8. The expression levels of genes encoding indole-3-acetic acid-amido synthetase (GH3), auxin-responsive Aux/IAA family members (IAA1, IAA20, IAA9), and auxin-responsive SAUR family members (SAUR11, SAUR33, SAUR25, SAUR55, SAUR31) were all significantly higher in a8 than b8 kernels. Two genes encoding auxin-repressed proteins (Aux rep) displayed significantly higher expression levels (> 10-fold) in a5 and a8 transcriptomes than in b5 and b8 transcriptomes. Of kernel development, the DEGs that encoded 6-phosphofructokinase (PFK), ESR6, ESR2, MRP-1, MEG-1, TCRR-1, ZmWRI1a, ZmAFL4, and wx1 had significantly higher expression levels in a8 than in b8 kernels. The transcript levels of invertase1 and colorless2 were significantly higher in ILb transcriptomes at all times than that in the respective ILa transcriptomes (Fig. 4, Additional file 6: Figure S5, Additional file 1: Table S1). These indicate that ILb kernels have higher sucrose metabolism as invertases hydrolyze sucrose into glucose and fructose.
We detected 53 genes that encode WRKY family members; 11 of these displayed significantly higher expression levels in a8 than in b8 kernels. One WRKY member (GRMZM2G029292) was specifically expressed in ILa kernels. PR10, PR1, and PR5 displayed significantly higher expression levels in a5 than in b5 kernels; PR5 expression level was > 10-fold higher in a5 kernels (276.72) than in b5 kernels (22.01) (Fig. 4). Other defense-related DEGs, including six peroxidase genes, four GST genes, and 18 CYP genes, all displayed significantly higher expression levels in a8 than in b8 kernels (Additional file 1: Table S1). These highly expressed DEGs involved in auxin, kernel development, or defense-response in a5 and a8 kernels may contribute to difference in early kernel development compared to b5 and b8 kernels and lead to smaller mature kernel size in ILa compared to ILb.
Some growth-related and defense-related genes showed dramatically abrupt changes in expression levels in ILa 8 DAP kernel. The expression levels of growth-related genes, encoding expansin, sucrose synthase1, et al., dramatically and abruptly decreased to almost zero in ILa 8 DAP kernels, and then recovered to similar levels as in ILb 10 DAP kernels. The expression levels of defense-related genes, encoding sulfur-rich/thionin-like protein, T22K18_16, hydrophobic protein RCI2B, senescence-associated protein DIN1, and natterin-4 dramatically, et al., abruptly increased in a8 kernels, and then decreased to similar levels as ILb kernels at 10 DAP (Additional file 7: Figure S6). These dramatically abrupt changes in expression levels suggest rapid global transcriptional reprogramming and an imbalance in growth and defense response in ILa kernels at 8 DAP, this might be associated with the simultaneous higher IAA and SA levels in ILa kernels at 8 DAP, compared to that of ILb kernels.
Dramatic transcriptional reprogramming occurred earlier in ILa than in ILb during early kernel development
ILa and ILb kernels displayed striking differences in their transcriptional profiles at 8 DAP, and they also showed their own specific transcription profiles from 5-DAP to 10-DAP. During kernel development from 5 to 8 DAP, the number of DEGs in ILa (1621 DEGs in the a8/a5 transcriptome pair) was more than twice the number of DEGs in ILb (716 DEGs in the b8/b5 transcriptome pair). The number of DEGs in the other two transcriptome pairs did not significantly differ; there were 1886 and 1937 DEGs in the a10/a8 and b10/b8 transcriptome pairs, respectively; and 2193 and 2207 DEGs in the a10/a5 and b10/b5 transcriptome pairs, respectively (Fig. 5a). These results indicate that dramatic transcriptional reprogramming in ILa kernels occurs during the transition from 5 to 8 DAP, whereas dramatic transcriptional reprogramming in ILb kernels occurs during the transition from 8 to 10 DAP. Within the ILa transcriptome pairs, a8/a5 shared more than half of its DEGs (835) with a10/a5, but only shared 485 of its DEGs with a10/a8, and a10/a5 shared 902 DEGs with a10/a8. The ILb transcriptome pairs had similar overlapping DEGs; b8/b5 shared 445 and 227 DEGs with b10/b5 and b10/b8, respectively, whereas b10/b5 shared 1425 DEGs with b10/8. During kernel development from 5 to 8 DAP, a8/a5 shared only 362 DEGs (22.33%) with b8/b5, with 1259 (77.67%) DEGs specifically owing to the a8/a5 transcriptome pair. From 8 to 10 DAP, transcriptome pair a10/a8 shared 695 (36.85%) DEGs with the b10/b8 transcriptome pair, each with almost two- thirds DEGs specific to its transcriptome pair, whereas a10/a5 shared more than half of its DEGs (1306) with b10/b5 from 5 to 10 DAP (Fig. 5b–d). A total of 111 TF-DEGs were differentially expressed in at least one of the three transcriptome pairs of ILa or ILb. We detected multiple DEGs encoding Hox and SBP TFs, whereas no WRKY members were differentially expressed in ILa or ILb at any developmental stage (Fig. 5e). These suggest that the greatest difference in transcriptional profiles between ILa and ILb is observed during the transition from 5 to 8 DAP. The dramatic transcriptional reprogramming induced a faster rate of developmental progression, which may begin at 5 DAP in ILa and 8 DAP in ILb, and this substantially different developmental processes ultimately contribute to morphological differences in ILa and ILb kernels.
The onset of the significantly enriched functional categories was earlier in ILa young kernel at 5 DAP but weaker at 10 DAP, compared to that in ILb young kernel
Functional analysis of DEGs from transcriptome pairs of ILa or ILb at different developmental stages showed that the a8/a5 transcriptome pair contained significantly more DEGs in all the top-enriched GO categories (biotic and abiotic stress responses) than the b8/b5 transcriptome pair. By contrast, and the enriched DEG number difference between a10/a8 and b10/b8 transcriptome pairs was not significant in most of these GO categories. Striking differences in the numbers of DEGs from ILa and ILb were observed in the GO category “response to auxin stimulus”, in which a8/a5 contained 44 DEGS, b8/b5 contained 15 DEGs; while a10/a8 contained 29 DEGs, and b10/b8 contained 130 DEGs. Similar trends were observed for “response to hormone stimulus”, “reproduction” and “post-embryonic development”. ILa transcriptome pairs (a10/a8 and a8/a5) contained more DEGs enriched in “chloroplast” than the respective ILb transcriptome pairs (b10/b8 and b8/b5) (Fig. 6a). Consistently, the a8/a5 transcriptome pair contained more DEGs that were enriched in all the top-enriched KEGG pathways than the b8/b5 transcriptome pair. The b10/b8 transcriptome pair contained more DEGs that were enriched in all these top-enriched KEGG pathways than the a10/a8 transcriptome pair. DEGs that were enriched in the functional categories of “starch and sucrose metabolism”, “plant hormone signal transduction”, and “phenylpropanoid biosynthesis”, were all significantly more in a8/a5 than that in b8/b5 transcriptome pairs; conversely, DEGs that were enriched in these functional categories were significantly more in b10/b8 than in a10/a8 transcriptome pairs. For the “starch and sucrose metabolism” category, a8/a5 contained 52 DEGs, b8/b5 contained 20 DEGs; a10/a8 contained 63 DEGs, while b10/b8 contained 79 DEGs; this suggests that starch and sucrose metabolic pathways are more active during 5 to 8 DAP and less active at later stages from 8 to 10 DAP in ILa kernels, compared to that in ILb kernels. For the “plant hormone signal transduction” category, a8/a5 contained 91 DEGs, b8/b5 contained 27 DEGs; a10/a8 contained 87 DEGs, while b10/b8 contained 93 DEGs (Fig. 6b). These results indicate that the onset of significantly enriched functional categories was earlier in ILa (at 5 DAP), whereas these functional categories were similarly or more significantly enriched in ILb at later stages (8 DAP).