Comparative transcriptome analysis identified candidate genes associated with kernel row number in maize

Plant materials and sample collection

Two maize inbred lines with distinct phenotypes in KRN, PHG35 and Dan598, were planted in the experimental field of Shandong Agricultural University (Taian, China). Dan598, an elite line and a parent of widely used commercial hybrids in China, originated from the Lüda Honggu mixed population and introduced hybrid varieties PN78599. The Expired Plant Variety Protection (ex-PVP) germplasm PHG35, derived from G3BD2 and H7FS6, was obtained from the North Central Regional Plant Introduction Station (USDA-ARS, Ames, IA, United States). The trial plants were planted in four-row plots at a density of 66,700 plants per hectare. Row spacing and plant spacing were about 0.6 and 0.25 m, respectively. The management during the entire growing period is the same as local field management practices. At least 30 plants were used for the measurement of KRN trait in each line. For RNA-Seq, young ear tissues from both inbred lines were collected at development stages V6 through V10 with three biological replications at each stage. At the V6 and V7 stages, 25–30 individuals were pooled per biological replicate; at the V8, V9, and V10 stages, 15–20 individuals were pooled for each replicate. The samples were immediately flash-frozen in liquid nitrogen and then stored at −80 °C for subsequent analysis.

Library construction, quality control, and read alignment

RNA isolation, cDNA library construction, and RNA sequencing were performed at Novogene Bioinformatics Technology (Beijing, China). Briefly, total RNA was extracted using a TRIzol reagent. The quality of extracted RNA was assessed using the NanoDrop ND-2000 spectrophotometer and an Agilent 2100 Bioanalyzer. After the RNA sample was qualified, mRNA was enriched using magnetic beads with Oligo(dT). Then, mRNA was fragmented and converted into single-stranded cDNA with random hexamers. Subsequently, the double-stranded cDNA was synthesized and purified with AMPure XP beads. The purified double-stranded cDNA underwent end repair, A-tailing, and adapter ligation, followed by size selection using AMPure XP beads. Finally, PCR amplification was performed, and the PCR products were purified using AMPure XP beads to obtain the final library. After library construction, a preliminary quantification was performed using Qubit 2.0, and the library was diluted. Then, the insert size of the library was detected using Agilent 2100, and accurate concentration was quantified using the Q-PCR method. All 30 cDNA libraries were sequenced on an Illumina NovaSeq 6,000 platform with a paired-end 150 bp sequencing strategy. The raw reads underwent initial filtering and quality control processes to generate clean reads using fastp (Chen et al., 2018). The filtered high-quality clean reads of each sample were aligned to the B73 reference genome (https://download.maizegdb.org/Zm-B73-REFERENCE-GRAMENE-4.0) using HISAT2 version 2.2.1 (Kim et al., 2019). Alignments in the SAM format were converted into BAM format using SAMtools version 1.6 (Li et al., 2009). Transcript abundance was estimated as Fragments per kilobase of transcript per million mapped reads (FPKM) and read counts using StringTie version 2.2.1 (Pertea et al., 2016). The statistical power of this experimental design was calculated using RNASeqPower (Hart et al., 2013).

Global transcriptome analyses and differential gene expression analysis

To assess and explore the global transcriptomic profiles, Principal component analysis (PCA) and Pearson correlation among the three biological replicates were performed (Hou et al., 2024; Zheng et al., 2021). Pairwise differential gene expression analyses were conducted by the DESeq2 package (version 1.41.12) based on read counts (Love, Huber & Anders, 2014). Genes with foldchange ≥ 2 and p-adjusted < 0.05 were considered differentially expressed genes (DEGs). DEGs identified by comparing the two inbred lines at the same developmental stages were termed line-specific DEGs, and DEGs identified between two consecutive developmental stages in the same line were termed development-specific DEGs. The gene expression patterns of DEGs were clustered by the Mfuzz package of R (Kumar & Futschik, 2007). Gene Ontology (GO) classification was performed using the clusterProfiler package (Wu et al., 2021).

Prediction of candidate genes underlying known QTLs and QTNs

To further delineate candidate genes underlying QTLs and QTNs, we compiled KRN-related QTLs (Yan et al., 2006; Ma et al., 2007; Guo et al., 2011; Lu et al., 2011; Cai et al., 2014; Liu et al., 2015b; Yang et al., 2015; Zhang et al., 2017; Fei et al., 2022) and QTNs (Liu et al., 2015b; Zhang et al., 2017; Li et al., 2018; An et al., 2020; Liu et al., 2020; Zhang et al., 2020; Zhang et al., 2022) from previous studies. Given that different reference genomes were used in these studies, the BLAST search was performed for the markers of these QTLs and QTNs against the B73 RefGen_v4 genome in maizeGDB. QTLs that overlapped across different studies (n ≥ 2) were merged into larger QTL intervals and were defined as QTL hotspots in this study. Subsequently, we mapped the line-specific DEGs in Phase I to QTL hotspots and the 200 kb genomic region surrounding a significant QTN (100 kb each upstream and downstream), and the genes in these regions were regarded as KRN-related candidate genes for further analysis. Transcription factor databases were obtained from PlantPAN 4.0 (http://plantpan.itps.ncku.edu.tw/).

Expression validation by quantitative real-time PCR (qRT-PCR)

A random sample from each of the three biological replicates at each stage of two inbred lines was used for reverse transcription and qRT-PCR analysis. Total RNA was extracted using the TRIzol method, and cDNA was synthesized using the Vazyme Reverse Transcription Kit R223. qRT-PCR was performed on an Applied Biosystems PCR instrument VIIA7 with the Ultra SYBR Mixture (Cwbio, Taizhou, China), following the manufacturer’s protocols. The 2-ΔΔCt method was employed to evaluate the relative abundance of the candidate gene, with three biological replicates. Expression levels were normalized to those of the endogenous control gene, CULLIN. The primer sequences can be found in Table S1.

Transcriptome sequencing of young ears in the two maize inbred lines

In this study, the average KRN for inbred line PHG35 was 20.5 ± 1.3 and was 12.1 ± 1.1 for inbred line Dan598 (Fig. 1). To explore transcriptional differences between the two maize inbred lines that have distinct KRN phenotypes, young ear tissues were collected from both inbred lines at development stages V6 through V10 with three biological replications each stage. A total of 30 cDNA libraries were constructed for sequencing. With a sample size of three, the average statistical power of this experimental design among different groups was 0.848. A summary of the sequencing data is shown in Table S2. An average of 45.2 million raw reads were generated for each library, and a total of 1,344.3 M high-quality clean reads were obtained after filtering (ranging from 40.13 M to 50.91 M per library; Q30 was 92.8% to 94.2%). The clean reads of each sample were then mapped to the maize B73 RefGen_v4 reference genome. The average mapping ratio was 89.40% (ranging from 87.99% to 90.96%) and the uniquely mapped reads ranged from 73.15% to 81.91%.

Global comparison of transcriptome among the samples

Pearson correlation analysis and principal component analysis (PCA) were performed on the 30 ear samples across various ear developmental stages in PHG35 and Dan598. Correlation coefficients were higher within each inbred line compared with that for different developmental stages (Fig. 2A). The two inbred lines were clearly distinct in PCA (Fig. 2B), suggesting that the two lines with distinct KRN had notably different transcriptomic profiles. Furthermore, the three early developmental stages (V6 to V8) were distinguishable from the two later stages (V9 and V10), suggesting that different transcriptional pathways could be activated at the different ear developmental stages. However, the sample D_V10_1 exhibited separation in Fig. 2B. Consequently, to ensure adequate data quality and reproducibility among replicates, we elected to exclude the sample D_V10_1. The Spearman’s rank correlation coefficients (SCC) between D_V10_2 and D_V10_3 were 0.97 (Fig. S1B), revealing sufficient data quality and reproducibility between the replicates. The average FPKM value of the biological replicate was used to determine the expression level for each stage in both lines. We found that 21,940, 21,558, 21,386, 21,804, and 22,258 genes were expressed (FPKM ≥ 1) at the V6, V7, V8, V9, and V10 stages, respectively, in PHG35; similarly, 21,436, 21,500, 21,417, 21,999 and 22,492 genes were found expressed at the corresponding stages in Dan598 (Fig. 2C). On average, 53.1% of the detected genes had FPKM < 1, and were considered not expressed; 22.3% had low expression (1 ≤ FPKM < 10), 22.6% were moderately expressed (10 ≤ FPKM < 100) and 2.0% were highly expressed (FPKM ≥ 100), without apparent difference in gene expression abundance distribution among different samples (Fig. 2D). The overall transcriptome dynamics of maize ear development in the two maize inbred lines are depicted in Fig. 2E.

Exploration of line-specific and development-specific differentially expressed genes

To identify line-specific DEGs, we conducted pairwise differential gene expression analyses between the PHG35 and Dan598 samples at each developmental stage. There were 3,497 DEGs (1,514 up-regulated and 1,983 down-regulated) at V6, 3,775 DEGs (2,501 up-regulated and 1,724 down-regulated) at V7, 3,508 DEGs (1,811 up-regulated and 1,697 down-regulated) at V8, 6,689 DEGs (3,651 up-regulated and 3,038 down-regulated) at V9, and 7,908 DEGs (4,244 up-regulated and 3664 down-regulated) at V10 detected in Dan598 when compared with PHG35 (Fig. 3A and Data S2). Taking all developmental stages together, we identified 11,897 DEGs between the two inbred lines, demonstrating immense differences in their transcriptomic profiles. Among these, V6-V8 stages contained 5,850 line-specific DEGs and V9-V10 stages contained 9,803 line-specific DEGs (Fig. 3B). Furthermore, we detected DEGs that were specifically expressed in one or more stages (Fig. S2). There were 454, 418, 310, 1,263, and 2,240 DEGs that were uniquely expressed in the V6, V7, V8, V9, and V10 stages, respectively. A unique set of 251 DEGs was common to both the V6 and V7 stages, including fasciated ear3 (fea3). Similarly, 226 DEGs were uniquely expressed in both the V7 and V8 stages, such as terminal ear1 (te1), barren stalk1 (ba1), and barren stalk fastigiate1 (baf1). Additionally, 353 common DEGs were specifically expressed from V6 to V8, such as barren inflorescence2 (bif2). The specific expression of 998 common DEGs spanned from V6 to V10, with examples including unbranched2 (ub2) and barren inflorescence4 (bif4). The genes ba1, baf1, bif2, ub2, and bif4 are involved in the initiation of lateral organs and axillary meristems (Chuck et al., 2014; Gallavotti et al., 2011; Gallavotti et al., 2004; Galli et al., 2015; McSteen & Hake, 2001), indicating the significance of DEGs during the V6-V8 period for maize ears formation. Between V9 and V10, 2,543 common DEGs were specifically expressed, while a few genes controlling maize inflorescence architecture were identified, like Zm00001d013603 (KRN5b, Shen et al., 2024).

To detect changes in gene expression levels at different ear development stages, we compared consecutive stages in each inbred line to identify development-specific DEGs. In PHG35, the expression of 503 genes (DEGs) significantly changed from V6 to V7, 499 DEGs were found from V7 to V8, and the DEG numbers became 1,291 and 1,576 from V8 to V9 and from V9 to V10, respectively (Fig. 3C and Data S3). The pattern of DEG occurrence in Dan598 was seemingly different, with 580 DEGs between V6 and V7, 346 between V7 and V8, 6,394 between V8 and V9, and 115 between V9 and V10 (Fig. 3C and Data S3). These results showed that the number of development-specific DEGs in Dan598 was much larger than in PHG35, particularly at the period from V8 to V9.

Expression patterns of DEGs during the ear development in the two inbred lines

As line-specific DEGs may better explain the variation in KRN between the two inbred lines, we performed Mfuzz clustering analysis on the 11,897 line-specific DEGs and identified four clusters in both Dan598 and PHG35. These clusters displayed diverse gene expression patterns (Figs. 3D–3K). Cluster 1 comprised 2,995 and 2,650 DEGs in Dan598 and PHG35, respectively, with expression levels gradually increasing from V6 to V8, peaking at V8, and then decreasing at V9 (Figs. 3D and 3H). Cluster 2 contained 2,289 and 3,614 DEGs in Dan598 and PHG35, respectively. The expression patterns showed low expression from V6 to V9 and then a significant increase at V10 (Figs. 3E and 3I). The expression patterns of the 3,618 DEGs in D_cluster 3 and the 2,836 DEGs in P_cluster 3 were lower from V6 to V8, showed an obvious increasing trend at V9, and then decreased at V10 (Figs. 3F and 3J). The expression level of the 2,995 DEGs in D_cluster 4 and the 2,797 DEGs in P_cluster 4 decreased gradually from V6 to V9 and slightly increased at V10 (Figs. 3G and 3K). Across the four clusters, V6 to V8 generally displayed different gene expression patterns compared to V9 and V10, suggesting an important developmental transition from V8 to V9, consistent with PCA clustering results. Previous studies have determined the vegetative growth stages by the leaf collar method (Abendroth et al., 2011; Ritchie, Hanway & Benson, 1986). Stages V6 to V8 begin 4 to 6 weeks after corn seedlings emerge from the soil. This phase is characterized by rapid growth and stem elongation, with the KRN being determined. Stages V9 and V10 mark the transition to a steady and rapid period of growth and dry matter accumulation. Therefore, we divided the five developmental stages into two distinct developmental phases: Phase I (V6 to V8) and Phase II (V9 and V10) in both PHG35 and Dan598.

We sought to examine whether there was a difference in GO terms between Phase I (V6 to V8) and Phase II (V9 and V10). Phase I contained 5,850 line-specific DEGs, and Gene ontology (GO) term enrichment analysis showed that these genes were significantly enriched in morphogenesis and differentiation processes, such as “regulation of plant organ formation”, “regulation of morphogenesis of a branching structure”, “plant organ formation”, “regulation of hormone levels” and “hormone transport” (Fig. 3L). Previous studies have indicated that the majority of genes affecting KRN are highly expressed in the early stage of phase I and are enriched in the aforementioned biological processes (Shen et al., 2023). The 9,803 DEGs between the two inbred lines in Phase II (from V9 to V10) were significantly enriched in growth-related terms like “cell cycle”, “DNA metabolic process”, “chromosome organization”, “nuclear division” and “meiotic cell cycle” (Fig. 3M). These findings indicate that Phase II represents a transition to rapid cell growth and cell proliferation, which corresponds to a key feature of ear development in this phase: the rapid expansion of ear coincides with the differentiation of floral organs (Abendroth et al., 2011; Ritchie, Hanway & Benson, 1986).

Mapping the line-specific DEGs in Phase I to KRN-related QTLs and QTN regions; prediction of candidate genes underlying known QTLs and QTNs

Considering the results of PCA clustering and GO term enrichment, we postulated that line-specific DEGs in Phase I are more likely involved in KRN determination in maize. Therefore, we mapped the 5,850 line-specific DEGs in Phase I to previously published KRN-related QTLs and QTN regions to further identify key regulatory genes for the KRN trait. A summary of KRN-related QTLs was shown in Data S4. LOD Score values of these QTLs range from 2.6 to 27.6, and PVE (%) range from 1.51 to 28.43. A total of 86 QTLs were identified across all ten chromosomes, and we depicted the 17 QTL hotspots (QTLs detected at least twice in previous studies) in Fig. 4, among which qKRN4.2, qKRN5.1, and qKRN2.1 had been detected 21, eight, and seven times, respectively; qKRN1.1, qKRN6.1, and qKRN10.1 were detected five times, and qKRN4.1, qKRN5.4, and qKRN7.1 were detected four times. A total of 1,999 line-specific DEGs in Phase I were included in the 17 QTL hotspots. Among these DEGs, several genes are known to affect ear development. For example, Terminal ear1 (te1, Zm00001d042445) was detected in qKRN3.1 and is considered a marker of ear initiation (White & Doebley, 1999); Zea floricaula/leafy1 (zfl1) and zfl2 were detected in qKRN10.1 and qKRN2.1, respectively. The zfl1 zfl2 double mutant led to a disruption of floral organ identity and patterning, and the activities of these two genes are significantly associated with KRN (Bomblies et al., 2003). Meanwhile, 347 QTNs significantly associated with KRN were identified from previously reported studies (Fig. 4 and Data S4). Within 100 kb upstream and downstream of these associated QTNs, 225 genes were line-specific DEGs in Phase I. In summary, there were 2,132 DEGs identified in the 17 QTL hotspots and the 347 QTN flanking regions combined, among which 92 DEGs were within QTLs where QTNs also exist, thus they were stated to be within the common regions of QTL+QTNs (Fig. 4 and Data S5).

The 92 DEGs within the common regions of QTL+QTNs may contain candidate genes that have a high impact on KRN. Except for qKRN2.2 and qKRN8.1, all the other QTL hotspots have a distribution of candidate genes (Fig. 4). Among them, qKRN5.2 and qKRN9.1 each contain one candidate gene, and qKRN4.2 contains the highest number of DEGs (17 genes). The expression trends of these genes were shown in boxes next to the gene names. In addition, nine transcription factors (TFs) were identified in the common regions of QTL+QTN. For example, Zm00001d002005, encoding an SBP-box TF, was located within the QTL qKRN2.1 with QTN snp55126 (Figs. 4 and 5B). Zm00001d052365 (AP2/EREBP 155) was mapped in QTL qKRN2.1 and QTN SYN3148, PZE-104109712 and SYN3156 regions (Figs. 4 and 5C). Zm00001d054012, encoding LRR-RLK, was located within QTL qKRN4.2 and QTN snp43631 region (Figs. 4 and 5E). Two hormone signaling pathway genes were identified: Zm00001d042267 (ARF10) was mapped to qKRN3.1 with two significant QTNs (SYN5822 and PZE-103097885), and Zm00001d025989 (Aux/IAA 43) was mapped to qKRN10.1 with two significantly QTNs (PZE-110080302 and PZE-110080486) (Figs. 4 and 5F).

Among the 2,132 DEGs that fall into QTLs or QTN regions, a total of 168 TFs were identified, including numerous development-related TFs, such as SBP-box, AP2/ERF, and MADS-box TFs (Fig. 5A). We identified four SBP-box genes, all of which were significantly up-regulated in Dan598 compared with PHG35 (Fig. 5B). Fifteen AP2/EREBP TFs were identified, with ten genes were down-regulated and five genes were up-regulated in Dan598 (Fig. 5C). Four MADS-box TFs were identified and all were down-regulated in Dan598 (Fig. 5D). We also identified nine leucine-rich repeat receptor-like protein kinase family proteins (LRR-RLKs), among which five were down-regulated in Dan598 (Fig. 5E). In addition, many genes involved in hormone homeostasis and signaling, such as auxin and cytokinin, were identified (Figs. 5F and 5G). The auxin-activating enzyme Zm00001d043701 (IAA-amino acid hydrolase ILR1-like 4) had higher gene expression in PHG35 at the V6 stage, while the auxin-inactivating gene Zm00001d010697 (IAA-amide synthetase GH3.6) had higher gene expression at the V8 stage in this line. Three auxin-responsive Aux/IAA family TFs (Aux/IAA 22, Aux/IAA 23 and Aux/IAA 37) were down-regulated in Dan598, whereas, Zm00001d034463 (Aux/IAA5) and Zm00001d025989 (Aux/IAA43) were up-regulated in Dan598 (Fig. 5F). CK-activating enzymes (LOG3, LOG5 and LOG8) and CK signaling component TYPE-A RESPONSE REGULATOR 16 (ARR16) were up-regulated in Dan598. Meanwhile, CK inactivation gene Zm00001d025885 (cytokinin oxidase11) was down-regulated in Dan598 compared with PHG35 (Fig. 5G).

qRT-PCR verification

To vindicate our transcriptome data, we selected eight genes to verify gene expression by quantitative real-time PCR (qRT-PCR; Fig. 6). By comparison, qRT-PCR results matched the RNA-seq data reasonably well for all these genes, indicating that our RNA-seq results are reliable and can be used for elucidating the regulatory networks for ear development and KRN determination.