Comparative transcriptomics analysis at the key stage of maize ear development dissect heterosis

Abstract

Important traits related to maize (Zea mays L.) grain yield, such as kernel row number, ear length, kernel number per row, are determined during the development of female inflorescence. There is a significant positive correlation between yield component and the activity of inflorescence meristem (IM). To find the key stage of heterosis in the development of the ear, immature ears (from the IM stage until the end of the floral meristem [FM] stage) of Yudan888 and its parent lines were sampled to assay phenotype and for comparative transcriptomics analysis. The immature ear length of Yudan888 at the IM stage fitted an additive (mid-parental) model, but it showed high parental dominance at the spikelet-pair meristem (SPM) stage. Comparative analysis of transcriptomes suggested significant differences between additive and nonadditive expression patterns for different developmental stages. The number of distinct maternal or paternal genes (DMP) (genes expressed only in one parental line and their hybrid but silenced in another line) was greater than ABF1 (genes expressed in both parental lines but silenced in hybrid) at each stage. Gene Ontology (GO) enrichment suggested that the cell redox homeostasis genes with overdominance expression patterns in hybrids have an important contribution to heterosis. According to our research, an ear length heterosis network was established. The discovery of the inflection point for ear length heterosis allows us for inferring that the transition state of IM to SPM may be the starting point of ear length heterosis. These findings improved the understanding of maize ear length heterosis.

Abbreviations
• ASE
• allele-specific expression
• DEGs
• differentially expressed genes
• FM
• floral meristem
• GO
• gene ontology
• GTP
• guanosine triphosphate
• IM
• inflorescence meristem
• qRT–PCR
• quantitative real-time–polymerase chain reaction
• SM
• spikelet meristem
• SPM
• spikelet-pair meristem

INTRODUCTION

Heterosis refers to the phenomenon that F1 individuals show stronger vigor, fertility, adaptability, stress resistance, and faster development than both parents (Birchler et al., 2010). In early 1908, George Shull first described the phenomenon of heterosis in maize (Zea mays L.) (Shull, 1908). Since then, geneticists have put forward several hypotheses, including dominance, overdominance and epistasis, to explain heterosis (Bruce, 1910; Cockerham, 1954; Davenport, 1908; Powers, 1945; Schnell & Cockerham, 1992). The dominance hypothesis holds that most dominant genes are conducive but that recessive genes are unfavorable; thus, dominant alleles in hybrids will typically dominate over recessive alleles, leading to heterosis (Bruce, 1910; Davenport, 1908). According to the overdominance hypothesis, interaction of alleles from both parents is the main cause of heterosis (East & Jones, 1920; Shull, 1908). The epistasis hypothesis emphasizes interaction between multiple loci (Cockerham, 1954; Schnell & Cockerham, 1992).

To improve the breeding efficiency of excellent hybrids, many studies on the genetic mechanism of heterosis have carried out. The quantitative trait loci (QTLs) linked to tillers/plant, grains/panicle, and seed weight were studied in a cross between Zhenshan 97 and Minghui 63, two elite rice (Oryza sativa L.) lines widely used in China. A large number of traits are related to digenic interactions, which provides evidence that epistasis is the main genetic basis of rice heterosis (Sibin et al., 1998). Using an immortalized F2, it shown that the dominant effect of heterotic loci at the single-locus level as well as AA (additive × additive) interactions plays significant roles in the inheritance of maize grain yield, ear length, ear row number, and 100-kernel weight (Tang et al., 2010). Several transcriptional sequencing studies detected mainly additive gene expression levels in hybrids (Guo et al., 2006; Thiemann et al., 2014; Zhou et al., 2019); others reported the prevalence of nonadditive gene expression (Ding et al., 2014; Fujimoto et al., 2012). Therefore, the generation of heterosis depends on species, traits, and genetic backgrounds, and no single hypothesis can fully explain the genetic mechanism of heterosis formation (Li et al., 2008). Some researchers believe that heterosis is not a consequence of higher levels of additive or nonadditive gene expression but of the change in transcription level between parents, with the transcriptional variation of specific genes affecting heterosis of specific traits (Stupar et al., 2008). A microRNA transcriptomic analysis showed approximately two-thirds of differentially expressed miRNAs in hybrids to be down-regulated. The target genes of the detected miRNAs, such as ARGONAUTE1 (AGO1), NUCLEAR FACTOR Y (NFY), SBP DOMAIN PROTEINS, and GROWTH REGULATING FACTORS (GRFs), also exhibited altered expression levels (Ding et al., 2012). This suggests that heterosis is a complicated trait affected by various types of regulatory components that act in hierarchies (Baranwal et al., 2012; Birchler et al., 2001, 2010; Veitia, 2010).

Important traits related to grain yield, such as kernel row number, ear length, and kernel number per row, are all determined during the development of the female inflorescence. An increase in inflorescence meristem (IM) diameter can provide more space for the differentiation of spikelet-pair meristem (SPM), which will lead to an increase in the number of rows per ear (Liu et al., 2015). The maintenance of IM activity can prolong the longitudinal development time of the ear, which makes it possible to accommodate more grains (Luo et al., 2021). Hence, there is a significant positive correlation between the yield component and differentiation activity of IM (Liu et al., 2015).

Yudan888, a new maize cultivar of our laboratory that has good grain quality, high disease resistance, and suitability for mechanical harvesting of kernels, is favored by growers. The mid-parent heterosis of Yudan888 is significant (Supplemental Figure S1). In this study, we continuously sampled the immature ears (from the IM stage until the end of the floral meristem [FM] stage) of Yudan888 and its parent lines (maternal inbred–15S717, paternal inbred–T4691). The ear length of Yudan888 at the IM stage fit an additive (mid-parental) model, but the SPM stage showed high parental dominance. Ear length shows an additive to dominant transition, which allows for inferring that differential expression of genes between the hybrid and parents at the SPM stage results in a heterotic phenotype and that the SPM stage is the best sampling time for studying ear heterosis. Therefore, we made comparative transcriptomics analysis at the key stage of maize ear development, which will provide new insights into the formation of maize ear-length heterosis.

MATERIALS AND METHODS Maize plants and field conditions

Yudan888 (15S717 × T4691) is a hybrid with wide adaptability, high yield, and good stability in the changeable environment of the Huang-Huai-Hai aera and is an ideal material for maize heterosis studies. In the summer of 2020, Yudan888 and its parental inbred Lines 15S717 and T4691 were grown at the Farm of Henan Agricultural University (Yuanyang, China). The planting density was 67,500 plants ha^–1, which were arranged in randomized blocks with three replications of 800 seedlings each.

Phenotyping

When the plants grew six (V6), seven (V7), eight (V8), and nine (V9) spreading leaves, the ear developed into the IM, SPM, spikelet meristem (SM), and FM stages, respectively. More than 100 ears of Yudan888 and its parental lines were collected at each stage for imaging under a stereomicroscope (Stemi 508, Carl Zeiss). When the ear is developing at the intersection of two phases, it is regarded as the representative sample of the previous stage (IM, SPM, SM stages). When a few bulges appear at the bottom of the IM, they indicate the transition from the IM to SPM stage. With this phenotype, ear length is determined as the IM stage length. The remarkable feature from SPM to SM is that spikelet lobes begin to appear in the middle of SPM and are divided into two pairs of spikelets, showing two rows of protrusions. When two pairs of spikelets appear at the bottom of the ear, the ear length is determined as the SPM stage ear length. When the outer glume of the spikelet appears, the ear length is determined as the SM stage ear length. To facilitate statistical analysis, we chose the stamen bulge appearing in the lower part of the immature ear as the key sign of the end of the FM stage (Cheng et al., 1983). The representative ear length was measured using ImageJ (Schindelin et al., 2015).

Core Ideas

• The spikelet-pair meristem period is the inflection point of heterosis of ear length.
• The distinct maternal or paternal genes may produce epistasis, leading to a nonadditive expression pattern.
• Nonadditive expression made the hybrid show stronger vigor than both parents.
• Cell redox homeostasis genes that maintain overdominant expression may have an important contribution to heterosis.
• An ear length heterosis network is described according to the study.

Developing ear samples

Inflorescence meristem tissue was dissected from the ear at four stages (IM, pre-SPM, post-SPM, SM stage) under a stereomicroscope (Stemi 305, Carl Zeiss) (Figure 1). Because an additive-to-dominance transition occurs during the SPM stage, this stage was further subdivided into pre-SPM (Figure 1b) and post-SPM stages (Figure 1c). When a few bulges appear at the bottom of the IM, the IM was cut from the ear as IM stage sample for transcriptomics analysis. When two pairs of spikelets appear at the bottom of the ear, the IM was cut from the ear as post-SPM stage sample for transcriptomics analysis. The phenotype of the ear between IM stage and post-SPM stage was classified as pre-SPM stage. When the outer glume of the spikelet appears, the IM was cut from the ear as SM stage sample for transcriptomics analysis. Three biological replicates were taken from developing ears of Yudan888 and its parents at each stage, with each biological replicate containing at least 30 IM ears. In total, 36 IM samples (4 developmental stages × 3 genotypes × 3 biological replicates) were collected for dissection of ear length heterosis. All samples were put into corresponding 1.5-ml Eppendorf tubes and stored at −80°C for further use.

FIGURE 1. Inflorescence meristem (IM) tissues at different inflorescence development stages were collected for assessing maize ear length heterosis. (a) IM stage. (b) Pre-spikelet-pair meristem (SPM) stage. (c) Post-SPM stage. (d) Spikelet meristem (SM) stage, Scale bars, 500 μm. The related meristem of each developmental stage is marked in the figure. The dotted line in the figure shows the IM cutting line under a stereomicroscope

RNA extraction, RNA-seq, and data analysis

RNA was extracted according to the TRIzol Reagent (Invitrogen/Ambion 15596-026) instructions. A total of 1 μg of RNA from each sample was used for the construction of RNA-seq libraries. Sequencing libraries were generated using NEBNextR Ultra TM Directional RNA Library Prep Kit for IlluminaR (NEB) The library quality was examined using an Agilent Bioanalyzer 2100 system (Agilent Technologies). The qualified libraries were sequenced using an Illumina HiSeq 2500 (HiSeq 2500, Illumina).

The adaptor and low-quality reads were filtered from the datasets by Trimmomatic (Bolger et al., 2014), and the clean reads were mapped to the B73 maize reference genome (http://ftp.ensemblgenomes.org/pub/release-40/plants/fasta/zea_mays/dna/) by HISAT2 (Kim et al., 2015). Only reads with a perfect match or one mismatch were further analyzed and annotated based on the reference genome. Differential expression analysis between sample pairs was performed using the DESeq2. Gene expression was divided into six categories (Figure 2) (Ding et al., 2012; Stupar et al., 2008). Gene Ontology (GO) enrichment analysis of the differentially expressed genes (DEGs) was implemented in the clusterProfiler R package (Yu et al., 2012).

FIGURE 2. Schematic diagram of hybrid gene expression categories. Lower than 10% of the expression level of the low parent is defined as under-dominant expression; expression between −10 and 10% of the low parent is defined as low-parent dominant expression; expression between −10 and 10% of the mid-parental value (MPV) is defined as additive expression; the expression difference between −10 and 10% of the high parent is defined as high-parent dominant expression; more than 10% of the expression level of the high parent is defined as overdominant expression; and other cases that do not conform to the above expression pattern are defined as “other”

Analysis of allele-specific expression

Allele-specific biases were evaluated using a binomial test with the null hypothesis that two parental alleles are uniformly expressed or modified in the hybrid. For allele-specific expression (ASE) analysis, single nucleotide polymorphism (SNP) calling was performed with GATK2 (McKenna et al., 2010). The GATK2 recognition criteria are as follows: there are no more than three consecutive single-base mismatches within the 35 bp range, and the SNP quality value is >2.0 after sequence depth normalization. The SNPs heterozygous in the offspring and homozygous in the parents were selected as candidates. Further analysis was performed when SNPs were identified with a significant allele-specific bias at a P value cut-off of .05. A gene was treated as an ASE gene only when candidate SNPs showed a significant difference (P < .05) in sequencing reads between the parents and offspring and the expression trend of all SNP sites in three biological repetitions was the same.

Validation of RNA-seq DEGs

The first-strand complementary DNA (cDNA) of the 36 RNA samples, which was used for RNA-seq, was synthesized with All-in-One First-Strand Synthesis MasterMix (KR0501-50 Kemix). Real-time fluorescence quantitative polymerase chain reaction (PCR) with 2× SYBR Green Master Mix (KS0601-500 Kemix) and a CFX96TM Real-Time System (Bio–Rad) was performed to quantify the expression levels of candidate genes. Three technical replicates were set in each plate for quantitative real-time–polymerase chain reaction (qRT–PCR). Maize ACTIN (Zm00001d010159) was used as the internal control. The 2-ΔΔCt method was used to evaluate relative gene expression levels (Livak & Schmittgen, 2001). The primers used for qRT- PCR are listed in Supplemental Table S1.

RESULTS Phenotypes of samples

Maize ears originate from meristems at the tip of lateral shoots (Cheng et al., 1983). When the plants had grown six spreading leaves (V6, 30 d after sowing for the hybrid and 33 d after sowing for both parental lines), the lateral shoot apical meristem (SAM) became the ear inflorescence meristem (IM) (Figure 3a). The IM length of Yudan888 was intermediate between that of the two parents (15S717, T4691) at the IM stage, adhering to an additive (mid-parental) model (Figure 3a and 3e). Next, a few bulges appeared at the bottom of the IM, which indicated the transition from the IM to SPM stage (Figure 3b). Yudan888 developed rapidly at this stage. Although there was no significant difference in the length of the immature ear between Yudan888 and T4691, the IM length of Yudan888 exceeded that of the two parents at the SPM stage and the difference was very significant (Supplemental Figure S2, Supplemental Table S2). When the plants develop to eight spreading leaves (V8), the ear develops into the SM stage. The remarkable feature from SPM to SM is that spikelet lobes begin to appear in the middle of SPM, showing two rows of protrusions (Figure 3c). At this stage, the ear length of Yudan888 exceeded that of the parent with the longer ear (T4691), and the difference was very significant (Figure 3h). That is, the ear length at the SM stage fit the overdominant model. The IM length of Yudan888 was consistently longer than both parents from the SPM to FM stage (Supplemental Figure S2). Thus, the immature ear length of Yudan888 showed an additive to overdominant transition during ear developmental processes. The change in IM length during ear developmental processes demonstrated a similar trend as the change in immature ear length (Supplemental Figure S2).

FIGURE 3. The meristem identities during inflorescence development of the three genotypes. (a–d) Micrographs of immature ears at the inflorescence meristem (IM) (a), spikelet-pair meristem (SPM) (b), spikelet meristem (SM) (c), and floral meristem (FM) (d) stages of the three genotypes. The yellow line represents the measured ear length, and the red line represents the measured IM length. Scale bars, 500 μm. (e–i) The immature ear length of Yudan888 and its parental lines at four inflorescence developmental stages. Single and double asterisks indicate significant differences in ear length at P [less than] .05 and P [less than] .01, respectively. “ns” indicates that the difference was not significant. M, maternal line (15S717); H, hybrid (Yudan888); P, paternal line (T4691); MPV, mid-parent value, which is a hypothetical value

RNA-seq analysis

Transcriptome sequencing of 36 libraries (4 developmental stages × 3 genotypes × 3 biological replicates) was performed. Approximately 43–52 million clean reads were obtained for each sample. The clean reads of each sample were mapped to the maize B73 reference genome, and the alignment efficiency ranged from 91.5 to 93.2%. The average correlation among the three biological replicates of each genotype at each stage ranged from 96.9 to 99.3% (Supplemental Table S3). The average fragments per kilobase of exon per million mapped fragments (FPKM) of the three replicates was taken as the expression level of the sample at each stage. In the four developmental stages, the number of genes shared by the three genotypes was 37,427, 37,277, 37,621, and 37,564. Screening was performed according to a certain standard (FPKM of hybrid≥1 or FPKM of mid-parent≥1), and the numbers after filtering were 21,779, 21,736, 22,079 and 22,168 for the IM stage, pre-SPM stage, post-SPM stage, and SM stage, respectively.

First, we analyzed the gene expression numbers of three genotypes at four developmental stages. As shown in previous studies (Wang et al., 2007; Xiong et al., 1999), five distinct types of differential gene expression patterns between hybrid and its parental lines appeared in this study, which included (a) UNF1: genes expressed only in hybrid; (b) DMP: genes expressed only in one parental line and their hybrid but silenced in another line; (c) ABF1: genes expressed in both parental lines but silenced in hybrid; (d) UNP: genes expressed in one parental line but silenced in another line and their hybrid; (e) MONO: genes expressed in both parental lines and hybrid (Figure 4). MONO type accounts for the largest proportion in each stage, reaching from 94.8 to 96.2%. The UNF1 and UNP, which are only specifically expressed in one genotype, were low at the four stages. Interesting, the number of DMP was greater than ABF1 at each stage (Figure 4a–d). Hence, the genes expressed in hybrids are more abundant than those of maternal or paternal lines.

FIGURE 4. Comparisons of gene expression types among the three genotypes. (a–d) Venn diagram comparisons among the expressed genes of the three genotypes at four developmental stages. (a) Inflorescence meristem (IM) stage; (b) pre-spikelet-pair meristem (SPM) stage; (c) post-SPM stage; (d) spikelet meristem (SM) stage. M, maternal line; P, paternal line; H, hybrid. (e–f) Multiple comparisons of co-expression between Yudan888 and (e) 15S717 or (f) T4691

We are interested in DMP type, which should belong to the allele-specific expression genes. Multiple comparisons of co-expression between Yudan888 and one parent (15S717 or T4691) at four stages were performed (Figure 4e and 4f). At the early stage of ear development, at least 115 genes were not expressed in T4691, and 261 genes were not expressed in 15S717, but they were all expressed in Yudan888. Gene annotation showed that these specific expression genes mainly include transcription factors and transposases. The percentage of transcription factors range from 10.53 to 17.06% (Supplemental Table S4). The DMP transcription factors between Yudan888 and one parent (15S717 or T4691) at IM stage are listed in Supplemental Tables S5–S6. These transcription factors can produce new interactions, which is also known as epistasis, and lead to changes in the quantity and quality of downstream gene expression and heterotic phenotypes of hybrids (Springer & Stupar, 2007).

We compared expression differences between the hybrid and parents at four developmental stages (Figure 5). Gene expression was divided into six categories: underdominance, low parental dominance, additive, high parental dominance, overdominance, other (see methods). The divergence in the total types of expressed genes in the hybrid was not significant at different developmental stages, ranging from 21,736 to 22,168. However, the additive/nonadditive expression gene numbers were significantly different (Figure 5). Additive expression was 51.8, 29.3, 55.4, and 21.8% at the IM stage, pre-SPM stage, post-SPM stage and SM stage, respectively. The overdominance (21.9%) and high parental dominance (15.2%) shown in the pre-SPM stage were higher than the same pattern in the other stages. During the SM stage, nonadditive expression increased significantly, mainly due to up-regulation of underdominance (32.6%). In the pre-SPM stage, overdominance is related to the heterotic phenotype of ear length.

FIGURE 5. Additive/nonadditive expression of the hybrid at four developmental stages. IM, inflorescence meristem; SM, spikelet meristem; SPM, spikelet-pair meristem

Analysis of allele-specific expression

A total of 6,457, 5,170, 6,718, and 6,583 ASEs were identified at the IM, pre-SPM, post-SPM and SM stages, respectively. In each period, the proportion of alleles favoring both parents was roughly equal (Figure 6a). A total of 9708 ASEs were identified in the four periods (Figure 6b). A total of 3,074 (31.4%) ASEs were common at each developmental stage. There are some specific ASEs in each stage, suggesting that ASE is stage specific. The bias of most ASEs was consistent at every stage of development, and only 5.7% (554/9708) of ASEs changed, suggesting a significant parent-of-origin effect for the action of parental alleles in hybrids.

FIGURE 6. Allelic bias of gene expression. (a) Percentage of alleles in four stages. Orange indicates allele-specific bias for the maternal line (15S717), green indicates allele-specific bias for the paternal line (T4691). (b) Selective expression analysis of alleles at four stages. IM, inflorescence meristem; SM, spikelet meristem; SPM, spikelet-pair meristem

Gene ontology enrichment of overdominant and underdominant genes in hybrid

To understand ear length heterosis, Gene Ontology (GO)-enrichment analysis was performed, mainly on overdominant and underdominant genes (Figure 7, Supplemental Tables S7 and S8).

FIGURE 7. Top 10 of Gene Ontology (GO) terms enrichment at four stage. Gene Ontology enrichment of overdominant genes at inflorescence meristem (IM) (a), pre-spikelet-pair meristem (SPM) (c), post-SPM (e), and spikelet meristem (SM) (g) stages; GO enrichment of underdominant genes at IM (b), pre-SPM (d), post-SPM (f), and SM (h) stages

The IM stage is the initial stage of floral organ development. Morphological studies show that the length of hybrid IM performance fits an additive (mid-parental) model (Figure 3a). Thus, transcriptome analysis at the IM stage is of great significance for understanding the formation of heterosis. Gene Ontology term enrichment in overdominant genes of the hybrid suggests that biological processes related to multicellular organismal development, response to water, and response to auxin become activated (Figure 7a). According to the highly active multicellular organismal development of the hybrid, six of the YAB3-, YAB4-, YAB5-, and DROOPING LEAF (DL)-related genes showed overdominant expression at the IM stage. Eight primary auxin response genes encoding small auxin up RNA proteins (SAUR25, SAUR32 SAUR36, SAUR71) and an indole-3-acetic acid-induced protein ARG7 were highly up-regulated in the hybrid. In contrast, small guanosine triphosphate (GTPase)-mediated signal transduction, intracellular signal transduction, and DNA topological change were enriched processes expressed in underdominance (Figure 7b).

At the pre-SPM stage, gene transcription was very active, with 4,769 genes showing overdominant expression in hybrids (Figure 5). The GO enrichment showed that distinct biological processes are activated in hybrids compared with their parental lines (Figure 7c). Forty genes belonging to the small GTPase-mediated signal transduction process, which are expressed in underdominance at the IM stage, showed overdominance expression at the pre-SPM stage. Due to the vigorous metabolism and active biochemical reaction of hybrids, 31 genes that maintain cell redox homeostasis were enriched as overdominant in the hybrid, including GLUTAREDOXINS (GRXS2, GRXS4, GRXS5, GRXC3, GRXC9), ADENOSINE 5′-PHOSPHOSULFATE REDUCTASE (APR1), PHOSDUCIN-LIKE PROTEIN (PLP3B), THIOREDOXIN (TRXH2, TXNL1, TRXX, TRXM1), and NUCLEOREDOXINS (NRX1). The nucleoside metabolic process was also enriched as overdominance. The underdominance ratio was lowest at the pre-SPM stage (Figure 5). The GO-term enrichment showed that the guanosine tetraphosphate metabolic process was not activated in hybrids compared with the low parent (Figure 7d).

Compared with the pre-SPM stage, the number of overdominant expressed genes decreased, but expression of underdominance increased in the post-SPM stage (Figure 5). It is suggested that protein synthesis, folding, modification, and cell growth are enhanced at this stage but that gene transcription is relatively decreased for hybrid. The biological processes of DNA replication, polysaccharide catabolism, and nucleoside metabolism showed overdominant transcription (Figure 7e). In contrast, the expression levels of genes related to protein ubiquitination, signal transduction, and guanosine tetraphosphate metabolism in hybrid offspring were lower than those in the low parent (Figure 7f). The ubiquitination-related proteases were lower in hybrids, suggesting fewer misfolded proteins or that the duration of protein activity was longer.

In the SM stage, the ear developed rapidly, and the morphological changes were great (Figure 3c). Nonadditive expression accounted for 78.19%, which was highest among the four stages. The ratios of overdominance and underdominance were 14.43 and 32.64%, respectively. The processes of small GTPase-mediated signal transduction, cell redox homeostasis, and translational initiation showed overdominance (Figure 7g). The processes of signal transduction, protein deubiquitination, and phosphatidylinositol metabolic process displayed characteristic lower expression in the hybrid compared with the low parent (Figure 7h).

Potential genes contributing to ear length heterosis

Nonadditively expressed genes may be sources of heterosis (Auger et al., 2005). To identify potential genes that contribute to ear length heterosis, multiple comparisons of overdominant and underdominant genes were performed in hybrids during four stages.

Twenty-six genes showed overdominance in all four stages (Figure 8a and 8c). Three genes with sustained overdominant expression were selected as the most candidate genes, namely, RPL11 (Zm00001d010016), GA20ox1 (Zm00001d013725), and EIN3 (Zm00001d022530). RPL11, an important component of the large subunit of ribosomes, is responsible for recruitment of initiation, elongation, and release factors (Mazzoleni et al., 2015). RPL11 induces activation of p53 to ensure cell growth (Zheng et al., 2015). GA20ox1, a key enzyme that catalyses the penultimate reaction of gibberellin biosynthesis, promotes floral organ growth and increases the number of grains (Plackett et al., 2012; Wu et al., 2016). At the primary growth stage of rice, OsGA20ox1 is a major gene that promotes plant height and leaf sheath growth (Oikawa et al., 2004; Yano et al., 2012). EIN3, a transcription factor that binds to the DNA sites of ethylene-responsive genes, controls a multitude of downstream transcriptional cascades to accelerate plant growth (Zhao et al., 2021).

FIGURE 8. Identification of potential genes contributing to ear length heterosis. (a, b) Venn diagram comparisons among overdominant (a) and underdominant (b) genes in the hybrid at the four stages. (c, d) The expression profiles of the common genes in a and b. Different columns in c and d represent the samples of three genotypes at the four stages, and different rows represent potential genes contributing to ear length heterosis. Color represents the expression level log10 of genes in the samples. IM, inflorescence meristem; SM, spikelet meristem; SPM, spikelet-pair meristem

Fifty-one genes were enriched in four developmental stages and showed underdominance (Figure 8b and 8d). MSSP2 (Zm00001d016274), PAPD5 (Zm00001d036629), and SYNDECAN-4 (Zm00001d012156) were selected as candidate genes related to ear length heterosis. MSSP2, a myc single-strand binding protein, contains RNP consensus motifs and is thought to function as a cell cycle regulator (Niki et al., 2001). PAPD5 is a noncanonical polymerase that can reduce telomerase activity and inhibit cell replication (Nagpal et al., 2020). Knockdown of SYN4 promotes cell proliferation during zebrafish embryogenesis (Luo et al., 2016).

Validation of differentially expressed genes by quantitative polymerase chain reaction

Real-time quantitative PCR (RT-qPCR) were performed to verify the accuracy of RNA-seq. Six candidate genes were selected, and their relative expression levels were quantified. The relative expression of all six genes were consistent with the RNA-seq FPKM data (Supplemental Figures S3–S6), which showed the reliability of our RNA-seq results.

DISCUSSION Heterosis is trait and genotype specific

Compared with the parents, both the time from sowing to emergence and the growth vigor at the seedling stage of the hybrid showed strong heterosis. Additionally, the SAM became the ear IM of the hybrid 2 d earlier than in the parents. Surprisingly, the hybrids of the IM stage were not longer than the high parent, which was equal to the mid-parent value and fit an additive model. The ear length then gradually caught up with the high parent at post-SPM stage (Figure 3). For both SM and FM stage, the length of the hybrid exceeded that of the higher parent, and the difference was significant. The hybrids showed additive and overdominance in sequence. The SPM is the key transmission point of ear length heterosis (Supplemental Table S2). A certain time point of IM-SPM may be the starting point of ear length heterosis, which is the best sampling time for studying ear heterosis.

Epistasis inducing nonadditive expression

In recent years, a great deal of work has reported on heterosis at the transcriptional level (Auger et al., 2005; Dooner et al., 1991; Goff & Zhang, 2013; Guo et al., 2006; Guo et al., 2014; Lu et al., 2003; Sibin et al., 1998; Tang et al., 2010; Thiemann et al., 2014; Yu et al., 1997; Zhou et al., 2019). Some studies show that additive expression is the main contribution to heterosis (Guo et al., 2006; Zhou et al., 2019), and some show that nonadditive expression is more important (Auger et al., 2005), which includes dominance (Guo et al., 2014), overdominance (Goff & Zhang, 2013; Lu et al., 2003) and epistasis (Dooner et al., 1991; Sibin et al., 1998; Tang et al., 2010; Yu et al., 1997). According to our ear length phenotypic data, we performed RNA-seq of the IM tissue at four key developmental stages (IM, pre-SPM, post-SPM, and SM stages). Transcriptome data showed gene expression to be a dynamic process during maize ear development (Figure 5). At the IM stage, half of the genes were transcribed additively, which corresponded to the mid-parent phenotype of hybrids. When the immature ear developed into the pre-SPM stage, 70.67% of the genes were nonadditively expressed. At this stage, overdominance was prominent, accounting for 21.94% of the total expressed genes. What causes the nonadditive expression of hybrid genes? Maize is rich in allelic variation due to mutation and transposition. An Affymetrix microarray study suggested that ∼2.5% of the genes present on the array were detected in only one of two transcriptomes, either B73 or Mo17 (Stupar & Springer, 2006). The DMP types showed the same results (1.3−2.6%) (Figure 4a–d). The DMP types could be the result of cis-acting factors of the parental gene (Springer & Stupar, 2007). The specific co-expression between the hybrid and one parent contains some transcription factors (Supplemental Tables S2–S3). These transcription factors may produce new interactions (epistasis) and combine to activate or suppress expression of target alleles in the F₁ hybrids, leading to a nonadditive expression pattern. Maize B and Pl genes encode transcription factors, and their interaction can enhance expression of A1, A2, and Bz1, which control anthocyanin production (Dooner et al., 1991). An inbred line with a nonfunctional b or pl allele displays a green phenotype owing to low or absent expression of genes A1, A2, and Bz1, but a hybrid with B/b Pl/pl alleles has high expression levels for genes A1, A2, and Bz1 and a red phenotype (Springer & Stupar, 2007). Therefore, we speculate that interactions of the specific co-expressed proteins at the IM stage lead to overdominant expression at the pre-SPM stage in the hybrid. Therefore, ear length fits a high-parent dominance model at the post-SPM stage.

Relationship between biological process and heterosis

Immature ears develop rapidly, and the proportion of nonadditive expression in hybrids at each period is relatively large, such as the pre-SPM and SM stages. Enrichment of biological processes showed a great difference between the IM stage and the next three stages of hybrids in terms of both overdominance and underdominance. The YABBY family, which is related to multicellular organismal development, is actively expressed in hybrids at the IM stage. YABBY genes encode seed plant-specific transcription factors that play pivotal roles in the modulation of morphogenesis and development (Zhang et al., 2020). The rice YABBY4 gene regulates plant growth and development by modulating the gibberellin pathway (Yang et al., 2016). Maize YABBY genes DROOPING LEAF1 and DROOPING LEAF2 regulate ear floret development (Strable, & Vollbrecht, 2019). SMALL AUXIN UP RNA genes (SAURs) are early auxin-responsive genes that also showed overdominance at the IM stage. As one of the converging points of auxin and gibberellin signals, SAURs regulate cell elongation and transcription of calcium and iron homeostasis genes in dicotyledonous plants (Agarwal et al., 2019; Qiu et al., 2020; Stamm & Kumar, 2013). Such overdominance of multicellular organismal development and the primary auxin response may be related to the generation of heterosis, which is upstream of the heterosis signalling pathway of ear length.

In pre-SPM as well as post-SPM and SM stages, the enriched biological processes were similar (Figure 7). Up-regulation or down-regulation of these biological processes together produce heterosis. Intriguingly, small GTPase-mediated signal transduction is underdominant at the IM stage and overdominant at all later stages. This large contrast suggests that biological processes may contribute to ear heterosis.

Here, we focus on the cell redox homeostasis that is expressed in the overdominant model. In recent years, an increasing number of studies have shown that reversible oxidative modification may be a new molecular mechanism in signal transduction regulation and become an important method in addition to phosphorylation, glycosylation, and ubiquitination modification (Haddad, 2004). The known intracellular systems responsible for disulfide reduction mainly include the thioredoxin (TRX) system and glutaredoxin (GRX) system, which were expressed as overdominance in the immature ears of hybrids in our study. The TRXs and GRXs are involved in multiple processes, such as protection of biomolecules from gamma irradiation (Yoon et al., 2013), modulation of redox-dependent signalling cascades (J, 2003), CO₂ assimilation (Cheng et al., 2014), cell differentiation (Zheng et al., 2020), correct folding of proteins (Yuen et al., 2013), and carbon and nitrogen metabolism (Baumann & Juttner, 2002). Maintaining cell redox homeostasis is beneficial to cell division and growth (Aziz et al., 2017; Castellano & Sablowski, 2008). The TRXs and GRXs sustained overdominant expression in the immature ear of the hybrid, possibly playing an important role in ear length heterosis.

Protein ubiquitination and deubiquitination mainly involve degradation of misfolded or oxidated inactive proteins (Vierstra, 2009). Overdominant expression of redox components in hybrids results in a few misfolded and inactivated proteins. Therefore, enzymes related to protein degradation are expressed in hybrids with underdominance. Overall, hybrids use more energy for growth and development (Goff & Zhang, 2013).

An ear length heterosis network

According to our phenotype and RNA-seq data, a mechanistic model of maize ear length heterosis was constructed (Figure 9). Fifty-two percent of genes in the hybrid were expressed additively in the IM stage of the immature ear, which may be the result of cis-acting genes. Interaction between parental chromosomes can lead to new expression patterns. The special epistasis from DMP types, which does not exist in parents, may activate or suppress expression of target alleles in the F₁ hybrids, leading to a nonadditive expression pattern. YABBY family genes and SAURs showed overdominance in the IM stage. The YABBY gene TONGARI-BOUSHI1 is involved in lateral organ development and maintenance of meristem organization in the rice spikelet (Tanaka et al., 2012). Mutation results in inflorescences comprising naked branches without spikelets (Tanaka et al., 2017). It has been suggested that the YABBY transcription factor is related to SPM differentiation. SAURs are early auxin-responsive genes. SAURs regulate cell elongation by integrating auxin and gibberellin signals (Stamm & Kumar, 2013). Therefore, we placed phytohormones in the signal shunt. Some studies have shown that circadian clock genes are also related to heterosis (Ni et al., 2009), and these signals affect the gene expression of metabolic pathways. The efficiency of DNA replication, protein translation, energy generation, and vesicle transportation of the hybrid is higher than that of the parental inbred lines. Enzymes related to protein degradation are expressed in hybrids with underdominance. Therefore, hybrids use more energy for growth and development. The reduction system plays an important role in cell metabolism. Its overdominant expression ensures the longevity and efficiency of various enzymes in hybrids, which is conducive to heterosis.

FIGURE 9. Mechanism of heterosis generation. The “+” in biologic process indicates up-regulation and “–” down-regulation

ACKNOWLEDGMENTS

The authors would like to thank Huili Yang, Yanan Lin and Guoqiang Xu for assisting with this article. This work was supported by a grant from the National Natural Science Foundation of China (31971961, 31871641, 32060460) and the Science and Technique Foundation of Henan Province of China (202102110036).

AUTHOR CONTRIBUTIONS

Liangfa Wang: investigation; methodology; Writing – original draft. Juan Li: investigation. methodology; Yuan Lin: Software; Visualization. Kuntai Dang: Investigation. Jiong Wan: Software; Validation. Shujun Meng: Formal analysis. Xiaoqian Qiu: Supervision. Qiyue Wang: writing—review & editing. Liqin Mu: Project administration. Dong Ding: Conceptualization; writing—review & editing. Hongbing Luo: Data curation. Jihua Tang: Conceptualization; Funding acquisition; Resources.

CONFLICT OF INTEREST

The authors declare that they have no conflict of interests.