Next Article in Journal
Seasonality, Composition, and Antioxidant Capacity of Limonene/δ-3-Carene/(E)-Caryophyllene Schinus terebinthifolia Essential Oil Chemotype from the Brazilian Amazon: A Chemometric Approach
Next Article in Special Issue
The Past, Present, and Future of Wheat Dwarf Virus Management—A Review
Previous Article in Journal
GABA Application Enhances Drought Stress Tolerance in Wheat Seedlings (Triticum aestivum L.)
Previous Article in Special Issue
Identification of Potential QTLs Related to Grain Size in Rice
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Plants
Volume 12
Issue 13
10.3390/plants12132496
plants-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Case Report

Genetic Analysis for the Flag Leaf Heterosis of a Super-Hybrid Rice WFYT025 Combination Using RNA-Seq

by
Qin Cheng
1,2,
Shiying Huang
1,2,
Lan Lin
1,2,
Qi Zhong
1,2,
Tao Huang
1,2,
Haohua He
1,2,* and
Jianmin Bian
1,2,*
1
Key Laboratory of Crop Physiology, Ecology and Genetic Breeding, Ministry of Education, Jiangxi Agricultural University, Nanchang 330045, China
2
College of Agronomy, Jiangxi Agricultural University, Nanchang 330045, China
*
Authors to whom correspondence should be addressed.
Plants 2023, 12(13), 2496; https://doi.org/10.3390/plants12132496
Submission received: 9 May 2023 / Revised: 6 June 2023 / Accepted: 27 June 2023 / Published: 29 June 2023
(This article belongs to the Special Issue Genetic Basis of Yield and Yield Stability in Major Crops)
Browse Figures
Review Reports Versions Notes

Abstract

:
The photosynthetic capacity of flag leaf plays a key role in grain yield in rice. Nevertheless, there are few studies on the heterosis of the rice flag leaf. Therefore, this study focuses on investigating the genetic basis of heterosis for flag leaf in the indica super hybrid rice combination WFYT025 in China using a high-throughput next-generation RNA-seq strategy. We analyzed the gene expression of flag leaf in different environments and different time periods between WFYT025 and its female parent. After obtaining the gene expression profile of the flag leaf, we further investigated the gene regulatory network. Weighted gene expression network analysis (WGCNA) was used to identify the co-expressed gene sets, and a total of 5000 highly expressed genes were divided into 24 co-expression groups. In CHT025, we found 13 WRKY family transcription factors in SDGhps under the environment of early rice and 16 WRKY family genes in SDGhps of under the environment of middle rice. We found nine identical transcription factors in the two stages. Except for five reported TFs, the other four TFs might play an important role in heterosis for grain number and photosynthesis. Transcription factors such as WRKY3, WRKY68, and WRKY77 were found in both environments. To eliminate the influence of the environment, we examined the metabolic pathway with the same SDGhp (SSDGhp) in two environments. There were 312 SSDGhps in total. These SSDGhps mainly focused on the phosphorus metallic process, phosphorylation, plasma membrane, etc. These results provide resources for studying heterosis during super hybrid rice flag leaf development.

1. Introduction

Heterosis is a phenomenon in which hybrids are superior to their parental lines in economic traits, which is an important strategy for crop breeding. A significant increase in rice grain yield in recent decades has been attributed to the use of heterosis in hybrids. At the genomic level, three hypotheses, namely dominant complementarity, over-dominance [1] and epistasis [2], are often proposed to explain the mechanism of heterosis. Although heterosis is widely used to enhance crop productivity worldwide, its molecular basis remains enigmatic. With the development of molecular marker technology and QTL mapping, geneticists used population genetic analysis in combination with molecular markers for QTL mapping of heterosis. Nevertheless, the partitioning of dominant populations and genetic diversity cannot fundamentally explain the heterosis that occurs between the offspring hybrids and their parents [3].
To uncover heterosis in hybrid rice, researchers used omics methods to find a number of genetic loci related to heterosis [4]. High-throughput RNA sequencing (RNA-seq) for in-depth transcript analysis contributes significantly to a better understanding of heterosis in rice [5]. Researchers analyzed the root transcriptomes of super-hybrid rice cultivar Xieyou 9308 and its parents at tillering and heading stages. They found that transcripts were differentially expressed between WFYT025 and its parents (DGHP) and identified a group of potential candidate transcripts [6]. Scientists used RNA-Seq technology to analyze the transcriptome of two rice hybrids, i.e., Ajay (based on wild-abortive (WA)-cytoplasm) and Rajalaxmi (based on Kalinga-cytoplasm), and their respective parents at the panicle formation (PI) and grain filling (GF) stages. The results showed that a number of important transcription factors (TFs) are encoded by the identified distinct expression genes (DEGs) [7]. E et al. used RNA-seq data from two rice genotypes and their reciprocal hybrids and identified genes that are differentially expressed or allele-specifically expressed in the hybrids and their parents. They identified 45 ASE (allele-specific expression) transcription factors belonging to 17 families [8]. Shao et al. performed a genome-wide analysis of ASE by comparing the reading ratios of SNPs of parental alleles in an elite rice hybrid (Shanyou 63) and its parents [Zhenshan 97 (ZS97) and Minghui 63 (MH63)] using RNA-seq data from seedling shoots, flag leaves, and young panicles of plants grown under four environmental conditions [9]. The results suggest that consistent ASEGs can cause partial to complete dominance effects on the traits they regulate, while directionally shifting ASEGs can cause overdominance.
Increasing food production is a long-term goal of crop breeding to meet the demands of global food security [10]. In rice, the uppermost three leaves, especially the flag leaf, are the main source of carbohydrates supply for grain development [11]. It is reported that the flag leaf produces more than 50% of the carbohydrates for the grains. Leaf size and angle, as the main components of plant architecture, affect photosynthesis in the canopy and crop productivity [12,13]. Rice leaf expansion is largely determined by longitudinal cell division, cell elongation, and cell expansion. Some genes have been shown to be associated with rice leaf phenotype, such as NARROW LEAF 1 (NAL1) [14], NAL2 and NAL3 [15], NAL7 [16], NARROW AND ROLLED LEAF 1 (NRL1) [17], and ABNORMAL VASCULAR BUNDLES (AVB) [18]. However, there are few studies on the heterosis of the rice flag leaf. Therefore, it is necessary to investigate the heterosis of the rice flag leaf.
In this study, we focused on heterosis in the rice cultivar WFYT025, which is widely grown in China. The length, width, and leaf area of the flag leaf showed great differences between WFYT025 and its female parent. Therefore, we sampled the flag leaf of both WFYT025 and its two parents for high-throughput transcriptome sequencing to investigate genes related to photosynthesis or transpiration and seed development. Gene regulatory network analysis revealed several large sub-networks representing interactions between genes with similar expression profiles, which were termed co-expression modules. Uncovering the function of these transcripts may provide useful information for understanding the molecular mechanism of heterosis. Several major sub-networks were found to represent interactions between genes with similar expression profiles via gene regulatory network analysis. Co-expressed gene sets were identified by weighted gene co-expression network analysis (WGCNA), and a total of 5000 highly expressed genes were divided into 24 co-expression groups. We found nine identical transcription factors in both phases. Except for five reported TFs, the other four TFs might play important roles in grain number and photosynthetic heterosis.

2. Results

2.1. Phenotypic Analysis of WFYT025 and Its Parents

In this study, we examined flag leaf length, width, leaf area, and 1000 spikelet weight or 1000-grain weight (TGW) of WFYT025 and its two parents at 1 day post-anthesis (DPA) and 10 DPA under different environmental conditions (Table 1). In the early rice environment, WFYT025 exhibited greater leaf length, leaf width, leaf area, and Spikelets weight at 1 DPA compared with WFB and CHT025. However, the leaf length of WFYT025 was longer than that of CHT025 but shorter than that of WFB, and the TGW of WFYT025 was higher than that of CHT025 but lower than that of WFB at 10 DPA. In the middle rice environment, the photosynthetic rate (μm/m2/s) of WFB was 12.63, CHT025 32.08, and WFYT025 37.95 at 1 DPA, WFB 30.53, CHT025 36.50 and WFYT025 30.35 at 10 DPA. The transpiration rate (μm/m2/s) of WFB was 4.55, CHT025 11.7, and WFYT025 12.33 at 1 DPA, WFB 10.97, CHT025 10.03, and WFYT025 5.37 at 10 DPA. Compared to WFB and CHT025, WFYT025 had greater leaf length, leaf width, and leaf area at 1 DPA and 10 DPA.
Heterosis of the middle parent (MPH) and Heterosis of the higher parent (HPH) were estimated for flag leaf Heterosis (Table 1). MPH showed a negative effect on leaf length around early rice and TWK around middle rice at 10 DPA. Except for MPH mentioned above, MPH values for all traits varied from 5.3 to 43.4%. Most of the phenotype parameters of WFYT025 show MPH, although a few traits have high parental Heterosis HPH.

2.2. RNA Sequencing of WFYT025 and Its Parents

A total of 1725.99 million reads were generated, with an average of 47.94 million reads per sample. The total number of valid reads is 1559.21 million, and the average valid read per sample is 43.31 million. In each case, the Q20 fraction was over 99%, and the Q30 base percentage was over 95% (Supplementary Table S1). Therefore, the quality of the data was very high and met the requirements for further analysis.
To confirm the accuracy and reproducibility of the RNA-Seq results, 12 genes were selected for qRT-PCR with the specific primers (Supplementary Table S2). The validation results for the 12 genes are shown in Supplementary Figure S1. The qRT-PCR results were all in agreement with the RNA-Seq data. Thus, our transcriptome sequencing results were credible.

2.3. Analysis of Differentially Expressed Genes

To analyze gene expression at different stages under different environmental conditions, gene expression was mainly measured by FPKM (Fragments per Kilobase of exon model Per Million mapped reads) or RPKM (reads Per Kilobase of exon model Per Million mapped reads) (Figure 1A). In the leaves of CHT025 under the early rice environment, the number of DEGs (p < 0.05) at 1 DPA and 10 DPA totaled 3072, of which 1937 were up-regulated, and 1135 were down-regulated. At the same time, the number of DEGs in WFB leaves at 1 DPA and 10 DPA were 5039 in total, among which 3244 were up-regulated, and 1795 were downregulated. In WFYT025 leaves, the number of DEGs at 1 DPA and 10 DPA was 3757, of which 1937 were up-regulated, and 1135 were down-regulated (Figure 1B). In the environment of middle rice, the number of genes up-regulated in CHT025, WFB, and WFYT025 was 892, 1273, and 819, respectively, whereas the number of genes down-regulated in CHT025, WFB, and WFYT025 was 616, 1934, and 2196, respectively (Figure 1C).

2.4. Gene Ontology (GO) Annotation and KEGG Pathway Enrichment Analysis

GO enrichment analysis of these DEGs identified several significant biological processes at different stages and under different environmental conditions. In the environment of early rice, “chloroplast stroma” and “oxidation-reduction” processes were significantly enriched in CHT025 (Figure 2A). Meanwhile, “closed chloroplast”, “ATP binding”, “cytosol”, and “mitochondrion” were significantly enriched in WFB (Figure 2B). At the same time, “the cytosol”, “chloroplast”, and “oxidation-reduction” processes were significantly enriched in WFYT025 (Supplementary Figure S2A). In the environment of the middle rice, “protein serine/threonine kinase activity”, “protein phosphorylation”, “extracellular region”, and “plasma membrane” were significantly enriched in CHT025 (Supplementary Figure S2B). “Closed protein serine/threonine kinase activity”, “ATP binding”, “plasma membrane”, and “defense response” were significant enriched in WFB (Supplementary Figure S2C). At the same time, “ATP binding”, “integral component of the membrane”, “plasma membrane” and “protein serine/threonine kinase activity” were significantly enriched in WFYT025(Supplementary Figure S2D). This suggests that the substantial differences in flag leaf between 1 DPA and 10 DPA may be related to photosynthetic efficiency.
Functional enrichment analysis was performed for all these DEGs during rice seed development. A total of 131 signaling pathways were identified in 1502 DEGs (including the 1 DPA and 10 DPA of CHT025 in the early rice environment). KEGG enrichment analysis showed that the t20 most enriched pathways included “ubiquinone”, “other terpenoid-quinone biosynthesis”, “glycerolipid metabolism”, “taurine”, “hypotaurine metabolism” and “tryptophan metabolism,” etc. (Figure 3A). Meanwhile, 2663 DEGs in WFB were classified into 137 metabolic pathways, and the 20 most enriched metabolic pathways mainly focused on concentrated “starch”, “sucrose metabolism”, “fatty acid metabolism”, “nicotinate”, “nicotinamide metabolism”, “fatty acid elongation,” etc. (Figure 3B). In total, 1911 DEGs in WFYT025 were classified into 133 pathways, and the 20 most enriched pathways mainly focused on concentrated “ribosome”, “purine metabolism”, “nitrogen metabolism” and “ribosome biogenesis in eukaryotes,” etc. (Supplementary Figure S3). In the environment of middle Rice, the pathways “plant–pathogen interaction”, “cutin”, “suberine”, “wax biosynthesis”, “brassinosteroid biosynthesis”, and “flavonoid biosynthesis” were enriched in CHT025 between 1 DPA and 10 DPA. In contrast, the pathways “plant-pathogen interaction”, “diterpenoid biosynthesis”, “phenylalanine metabolism”, and “flavonoid biosynthesis” were enriched in WFB. While “plant-pathogen interaction”, “amino sugar”, “nucleotide sugar metabolism”, “diterpenoid biosynthesis”, and “AGE-RAGE signaling” pathway in diabetic complications were enriched in WFYT025.

2.5. Identification of Gene Co-Expression Modules in the Flag Leaf of WFYT025

After obtaining the gene expression profile of the flag leaf, we further investigated the gene regulatory network. Analysis of the gene regulatory network revealed several large subnetworks representing interactions between genes with similar expression profiles, which were termed co-expression modules [19]. The co-expressed gene sets via WGCNA were identified, and a total of 5000 highly expressed genes were divided into 24 co-expression groups (Figure 4). We then analyzed the correlation between the gene expression profile of the co-expression module and the three phenotypic data (TGW; length; width) (Supplementary Figure S4). Since the goal is to find genes related to grain development, we focused on the co-expression module with a high correlation with TGW (the color of the module is Midnight blue, and the gene expression profile in this module is moderately negatively correlated with the thousand-grain weight phenotype). We found the overall gene expression profile of the co-expression module Midnight Blue, which had the highest correlation with TGW, and the co-expression module contained a total of 106 genes (Supplementary Figure S5). The GO enrichment analysis of these genes revealed several significant biological processes. Modules associated with TKW showed biological processes that are associated with “protein amino acid phosphorylation”, “protein modification process”, and so on (Figure 5). Overall, these modules represent the specific gene regulatory processes at each stage and are the indicators of operated programs of flag leaf development.

3. Discussion

In the plant life cycle, leaves play a significant role in grain yield, and flag leaf length has been recognized as an important factor that determines plant type for high-yield potential in rice [20]. However, most of the previous studies on rice yield focused on rice grains, and few focused on the effect of rice flag leaf on rice yield. In the current study, we investigated the relationship between transcriptional profiles and heterosis in super hybrid rice WFYT025 by RNA-Seq.

3.1. The Genetic Basis of Heterosis

A significant difference gene (SDG) was defined as the significant difference gene between 1 DPA and 10 DPA. Meanwhile, SDGs that existed in hybrid but not in parents were defined as SDGhp. We have been able to identify a number of SDGhps underlying flag leaf between hybrid WFYT025 and paternal line CHT025, confirming the suggestion that heterosis is a polygenic phenomenon [21]. Under the environment of early rice, among the SDGhps, 10.9% had a dominant effect, 41.81% had a partial dominant effect, 22.07% had an additive effect, and the remaining 25.22% had an over-dominant effect at 1 DPA. Meanwhile, at 10 DPA, among the SDGhps, 12.28% had a dominant effect, 42.04% had a partial dominant effect, 18.52% had an additive effect, and the remaining 27.16% had an over-dominant effect. Thus, partial dominance was the major contributor to the flag leaf heterosis of WFYT025 (Figure 6)
However, at 1 DPA under the environment of middle rice, among the SDGhps, 43.19% were partially dominant, 16.65% were dominant, 32.17% were over-dominant, and 7.99% had an additive effect. In contrast, at 10 DPA, 52.83% of the genes were super dominant, 26.54% partial dominant, 9.56% dominant, and 11.07% additive. These results might be due to high temperature, which was inconsistent with the results of other environments (Figure 7).

3.2. Transcription Factors May Underlie Heterosis

It is well known that the analysis of expression patterns of important transcription factors is helpful in understanding the heterosis in rice. In CHT025, we found 13 transcription factors of the WRKY family in SDGhps under the environment of early rice and 16 genes of the WRKY family in SDGhps of under the environment of middle rice. Transcription factors, including WRKY3, WRKY68, and WRKY77, were found in both environments. Recently, studies have shown that WRKY family genes are involved in the process of seed growth and matching in Arabidopsis and rice [22]. However, the expression level of WRKY68 in WFYT025 and its parents at 10 DPA was significantly higher than that of 1 DPA. The expression level of WRKY68 in WFYT025 was significantly higher than that of its parents. In the two stages, we found nine identical transcription factors. Except for five reported TFs, the other four TFs might play an important role in grain number and photosynthesis heterosis.

3.3. The Key Metabolic Pathways of Heterosis

In order to eliminate the impact of the environment, we studied the metabolic pathway with the same SDGhp (SSDGhp) in two environments. There were 312 SSDGhps in total. These SSDGhps were mainly concentrated in the phosphorus metallic process, phosphorylation, plasma membrane, etc. During senescence in annual crop plants, P (phosphorus) is mobilized from leaves and other vegetative tissues and translocated to seed development, which is a strong P sink during the reproductive growth phase [23,24,25]. Previous studies have shown that seed development is a strong sink for P, and the remobilization of P from vegetative tissues results in insufficient P for other processes necessary for continued growth, such as photosynthesis. In rice (Oryzasativa L.), P is predominantly mobilized from leaves during the later stages of grain filling [26,27]. This might be the reason why the filling speed of rice was faster from 1 DPA to 10 DPA after flowering than from 15 DPA to 21 DPA.

4. Material and Methods

4.1. Plant Materials and Growth Conditions

The hybrid WFYT025, along with its parental lines Changhui T025 (CHT025) and Wufeng B (WFB), were planted in the experimental field of Jiangxi Agricultural University, Nanchang (28°450′ N, 115°500′ E) in 2019. The sowing time of early rice was on 17 March, and the sowing time of middle rice was on 20 May. WFYT025 is a super-hybrid rice combination derived from the cross between female parent WFB and male parent CHT025. Each plot consisted of 40 rows, with each row having 10 plants at a 20 cm planting space. Crop management followed normal procedures for rice. The three lines were chosen in this study to measure phenotypic traits and conduct transcriptome analyses. We sampled the grains and flag leaf on the first day, the tenth day, the fifteenth day, and the twenty-first day after the anthesis of rice. By measuring the TGW of the grains, we finally selected the flag leaf at 1DPA and 10 DPA for sequencing. The flag leaf was collected and stored at −80 °C for RNA-Seq analysis, and each sample had at least three biological replications to minimize systematic errors.

4.2. RNA Extraction, Sequencing, De Novo Assembly and Statistical Methods

Total RNA was extracted from rice leaf using Trizol reagent (Invitrogen, Carlsbad, CA, USA) and purified using an RNeasy Plant Mini Kit (Qiagen, Valencia, CA, USA) according to the manufacturer’s instructions. The integrity of RNA was verified by RNase-free agarose gel electrophoresis, and the concentration was measured with a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). cDNA library was constructed using the mRNA fragments as templates and sequenced on the Illumina HiSeq™ 4000 platform (Illumina Inc., San Diego, CA, USA). Adapter sequences were removed from the raw reads. It follows that low-quality reads with 50% bases with quality scores of five or lower and/or over 10% bases unknown were removed to gain more reliable results. The raw data generated by sequencing were preprocessed by following Lianchuan Biological cutadapt procedure to filter out unqualified sequences to obtain valid data (clean data) before proceeding to the next step of the analysis. Raw data were cleaned through the following processing steps: (1) removing the reads with connectors (Adaptor); (2) removing the reads containing N (N indicates that base information cannot be determined) with a percentage greater than 5%; (3) removing low quality reads (the number of bases with quality value Q ≤ 10 accounts for more than 20% of the whole read); (4) counting raw sequencing volume, effective sequencing volume, Q20, Q30, GC content, and comprehensive analysis.

4.3. Pathway Enrichment Analysis

KEGG (Analysis of Kyoto Encyclopedia of Genes and Genomes) is the major public pathway-related database. Pathway enrichment analysis identified significantly enriched metabolic pathways or signal transduction pathways in DEGs (different expression genes) compared with the whole genome background. Differential enrichment analyses, including Gene Ontology and KEGG pathway enrichment analysis, enrichment analysis using the cluster profile R software package (edgeR).Significantly enriched pathways in DEGs compared to the genome background were defined by a hypergeometric test. The calculated p-value was gone through FDR Correction, taking FDR ≤ 0.05 as a threshold. Pathways meeting this condition were defined as significantly enriched pathways in DEGs.

4.4. GO Enrichment Analysis

Gene Ontology (GO) is an international standardized gene functional classification system that offers a dynamic-updated controlled vocabulary and a strictly defined concept to comprehensively describe the properties of genes and their products in any organism. GO has three ontologies: molecular function, cellular component, and biological process. The basic unit of GO is GO-term. Each GO term belongs to a type of ontology.
GO enrichment analysis provides all GO terms that are significantly enriched in DEGs compared to the genome background, filtering the DEGs that correspond to biological functions. GO enrichment analysis was performed using the OmicShare tools, a free online platform for data analysis (www.omicshare.com/tools accessed on 1 May 2023). Firstly all DEGs were mapped to GO terms in the Gene Ontology database (http://www.geneontology.org/ accessed on 1 May 2023), gene numbers were calculated for every term, and significantly enriched GO terms in DEGs compared to the genome background were defined by hypergeometric test. The calculated p-value was gone through FDR Correction, taking FDR ≤ 0.05 as a threshold. GO terms meeting this condition were defined as significantly enriched GO terms in DEGs. This analysis was able to recognize the main biological functions that DEGs exercise.

4.5. Quantitative Real-Time PCR (qRT-PCR) Validation

To validate the RNA-seq results, different expression patterns of several genes were confirmed by quantitative real-time RT-PCR (qRT-PCR). For qRT-PCR, 1 µg of total RNA was used to synthesized cDNA using PrimeScriptTM RT reagent Kit (Perfect Real Time) (Dalian, China, TaKaRa). The qRT-PCR was carried out using SYBR® Premix Ex Taq II (Tli RNaseH Plus; TAKARA BIO Inc., Shiga, Japan) and determined in LightCycler 480 (Roche, Basel, Switzerland) according to the manufacturer’s instructions. The qRT-PCR reactions were amplified for 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s, 55 °C for 30 s and 72 °C for 30 s. All reactions were performed with three independent biological replicates for each sample, and three technical replicates for each biological replicate were analyzed. The relative gene expression was calculated by the software of ABI7500 Real-Time PCR System using the 2−∆∆Ct method. The detection of the threshold cycle for each reaction was normalized against the expression level of the rice Actin1 gene with the primer sequences 5′-TGGCATCTCTCAGCACATT CC-3′ and 5′-TGCACAATGGATGGGTCAGA-3′.

4.6. The Mode of Inheritance Analysis

For statistical analysis, the analysis of variance (ANOVA) was usually by the model: y = u + (GA) + (GD) + (SR) + e, where y is the acquired gene expression, u is the overall mean, GA is the additive effect, GD is the dominant effect, SR is the replication effect, and e is the residual error (Chen L et al., 2018). Hp = [d]/[a], referred to as the dominance ratio or potency (where [a] and [d] represent GA and GD, respectively). Considering gene expression levels as quantitative traits, we adopted traditional quantitative genetic parameters, such as composite additive effect [a] and composite dominance effect [d], to estimate our expression profile. DGHP were classified according to the dominance ratio Hp (=[d]/[a]), based on 99.8% confidence intervals constructed for [d] − [a] ([d] > 0) and [d] + [a] ([d] < 0). According to the value of Hp (=[d]/[a]), we considered that these genes belonged to partial dominance (−0.8 < Hp ≤ −0.2 or 0.2 < Hp ≤ 0.8), over-dominance (Hp ≤ −1.2 or Hp > 1.2), dominance (−1.2 < Hp ≤ −0.8 or 0.8 < Hp ≤ 1.2), and additive effect (−0.2 < Hp ≤ 0.2) [28,29].

4.7. Normalization of Gene Expression Levels and Identification of Differentially Expressed Genes

Sequencing reads were mapped to the reference sequences. The expression level of each gene was measured by fragments per kilobase of exon model per Million mapped reads (FPKM). To determine the time-dependent transcriptional differences between early rice and middle rice, the differential expression genes (DEGs) at 1 DPA and 10 DPA were determined by comparing the expression levels. To correct for multiple testing, the false discovery rate (FDR) was calculated to adjust the threshold of the p-value [30]. The standard of different expressions between compare is a minimal 2-fold difference in expression (|log2 Ratio| ≥ 1) and an FDR ≤ 0.005 [31].

4.8. Weighted Gene Co-Expression Network Analysis (WGCNA)

The gene co-expression network analysis used the raw expression counts (RAW counts) of all Rice reference genes in the RNA-seq data of the 36 leaf groups obtained from the previous analysis. Firstly, the genes with low expression were eliminated. The elimination standard was that if the genes were not expressed in more than 80% of the samples, the genes with low expression were considered. Then, the median absolute deviation was selected from the genes with high expression(mad) 5000 genes of the maximum value for total express network analysis using WGCNA (https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/ accessed on 1 May 2023). In the process of WGCNA analysis, 20 soft power was selected. GO enrichment analysis through the agriGOv2 (http://systemsbiology.cau.edu.cn/agriGOv2/ accessed on 1 May 2023), its parameters as follows: Statistical Test Method as HyperGeometry, Multi_ Test Adjustment Method as Hochberg FDR, Significance level as “0.05” and a Minimum number of mapping entries as “5”.

5. Conclusions

The flag leaf of WFYT025 and its two parents were sampled for high-throughput transcriptome sequencing in identifying genes related to leaf photosynthesis or transpiration and the development of seeds. We selected the flag of hybrid rice WFYT025 and its parents for transcriptome sequencing, the 1 DPA and 10 DPA in the environment of early rice and middle rice, respectively. Several major sub-networks were found to represent interactions among genes with similar expression profiles via gene regulatory network analysis. The co-expressed gene sets via weighted gene co-expression network analysis (WGCNA) were identified, and a total of 5000 highly expressed genes were divided into 24 co-expression groups. In the two stages, we found nine identical transcription factors. Except for five reported TFs, the other four TFs may play an important role in grain number and photosynthesis heterosis.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/plants12132496/s1, Figure S1: Comparison of the log2 (FC) of 12 randomly selected transcripts using RNA-Seq and qRT-PCR. A–C: In the early rice environment, compare 1 DPA and 10 DPA about CHT025, WFB and WFYT025, respectively. D–E: In the middle rice environment, compare 1 DPA and 10 DPA about CHT025, WFB and WFYT025, respectively; Figure S2: Gene ontology (GO) enrichment analysis of DEGs. (A) GO enrichment analysis of DEGs at 1 DPA and 10 DPA under the environment of early rice in WFYT025. (B) GO enrichment analysis of DEGs at 1 DPA and 10 DPA under the environment of middle rice in CHT025. (C) GO enrichment analysis of DEGs at 1 DPA and 10 DPA under the environment of middle rice in WFB. (D) GO enrichment analysis of DEGs at 1 DPA and 10 DPA under the environment of middle rice in WFYT025; Figure S3: KEGG analysis of DEGs between WFYT025 in the environment of early rice. That showed the top 20 most represented categories and the number of transcripts predicted to belong to each category; Figure S4: Correlation between gene expression profile of co-expression module and three phenotypic data (TGW; Length; Width). The legend on the right shows the correlation coefficients, with positive correlation in red and negative correlation in blue. Since the aim is to find genes related to seed development, the focus is on the co-expression module with high correlation with TGW (module color is midnightblue, the gene expression profile in this module is moderately negatively correlated with the TGW phenotype); Figure S5: Is the overall gene expression profile of the co-expression module MidnightBlue, which has the highest correlation with TGW. There are 106 genes in this co-expression module. The graph is divided into three parts, the upper part is the expression heat map of the genes in the co-expression module. The middle part shows the relative expression of the eigengene in the module. The lower part shows the TGW phenotype values of the samples; Table S1: Summary of reads quality control; Table S2: Primers used in this study.

Author Contributions

Q.C. designed and performed experiments, analyzed data, and wrote the manuscript. S.H. participated in designing and performing experiment. L.L., H.H. and J.B. completed the manuscript with inputs in technical support, critical writing, and suggestions regarding the manuscript, Q.Z. and T.H. participated in performing experiments. All authors have read and agreed to the published version of the manuscript.

Funding

Jiangxi double thousand plan (jxsq2023201057); Key R&D Plan of Jiangxi Province (20224BBF62001; 20224BBF61030); Joint research on improved rice varieties of Jiangxi (JXNZWZY01).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are available upon request.

Acknowledgments

We thank the anonymous referees for their critical comments on this manuscript.

Conflicts of Interest

The authors declare that they have no conflict of interest.

References

  1. Birchler, J.A.; Yao, H.; Chudalayandi, S.; Vaiman, D.; Veitia, R.A. Heterosis. Plant Cell 2010, 22, 2105–2112. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  2. Kumar, A.; Sharma, V.; Jain, B.T.; Kaushik, P. Heterosis Breeding in Eggplant (Solanum melongena L.): Gains and Provocations. Plants 2020, 9, 403. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  3. Bian, J.M.; Jiang, L.; Liu, L.L.; Xiao, Y.H.; Wang, Z.Q.; Zhao, Z.G.; Zhai, H.Q. Identification of Chromosome Segments Associated with Heterosis for Yield in × Rice Hybrids. J. Rheumatol. 2010, 37, 1673–1679. [Google Scholar]
  4. Fu, D.; Xiao, M.; Hayward, A.; Jiang, G.; Zhu, L.; Zhou, Q.; Li, J.; Zhang, M. What is crop heterosis: New insights into an old topic. J. Appl. Genet. 2015, 56, 1–13. [Google Scholar] [CrossRef]
  5. Ge, X.; Chen, W.; Song, S.; Wang, W.; Hu, S.; Yu, J. Transcriptomic profiling of mature embryo from an elite super-hybrid rice LYP9 and its parental lines. BMC Plant Biol. 2008, 8, 114. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Zhai, R.; Feng, Y.; Wang, H.; Zhan, X.; Shen, X.; Wu, W.; Zhang, Y.; Chen, D.; Dai, G.; Yang, Z.; et al. Transcriptome analysis of rice root heterosis by RNA-Seq. BMC Genom. 2013, 14, 19. [Google Scholar] [CrossRef] [Green Version]
  7. Katara, J.L.; Verma, R.L.; Parida, M.; Ngangkham, U.; Molla, K.A.; Barbadikar, K.M.; Mukherjee, M.; Samantaray, S.; Ravi, N.R.; Singh, O.N.; et al. Differential Expression of Genes at Panicle Initiation and Grain Filling Stages Implied in Heterosis of Rice Hybrids. Int. J. Mol. Sci. 2020, 21, 1080. [Google Scholar] [CrossRef] [Green Version]
  8. Huang, S.; Zhang, Y.; Ge, L.; Wang, L. Genome-wide transcriptome profiles of rice hybrids and their parents. Int. J. Mol. Sci. 2014, 15, 20833–20845. [Google Scholar]
  9. Shao, L.; Xing, F.; Xu, C.H.; Zhang, Q.H.; Che, J.; Wang, X.M.; Song, J.M.; Li, X.H.; Xiao, J.H.; Chen, L.L.; et al. Patterns of genome-wide allele-specific expression in hybrid rice and the implications on the genetic basis of heterosis. Proc. Natl. Acad. Sci. USA 2019, 116, 5653–5658. [Google Scholar] [CrossRef] [Green Version]
  10. Huang, X.H.; Yang, S.H.; Gong, J.Y.; Zhao, Q.; Feng, Q.; Zhan, Q.L.; Zhao, Y.; Li, W.J.; Cheng, B.Y.; Xia, J.H.; et al. Genomic architecture of heterosis for yield traits in rice. Nature 2016, 537, 629–633. [Google Scholar] [CrossRef]
  11. Sanchez-Bragado, R.; Vicente, R.; Molero, G.; Serret, M.D.; Maydup, M.L.; Araus, J.L. New avenues for increasing yield and stability in C3 cereals: Exploring ear photosynthesis. Curr. Opin. Plant Biol. 2020, 56, 223–234. [Google Scholar] [CrossRef] [PubMed]
  12. Cui, K.H.; Peng, S.B.; Xing, Y.Z.; Yu, S.B.; Xu, C.G.; Zhang, Q. Molecular dissection of the genetic relationships of source, sink and transport tissue with yield traits in rice. Tag. Theor. Appl. Genet. 2003, 106, 649. [Google Scholar] [CrossRef] [PubMed]
  13. Wang, P.; Zhou, G.L.; Yu, H.H.; Yu, S.B. Fine mapping a major QTL for flag leaf size and yield-related traits in rice. Theor. Appl. Genet. 2011, 123, 1319–1330. [Google Scholar] [CrossRef] [PubMed]
  14. Qi, J.; Qian, Q.; Bu, Q.Y.; Li, S.Y.; Chen, Q.; Sun, J.Q.; Liang, W.X.; Zhou, Y.H.; Chu, C.C.; Li, X.G.; et al. Mutation of the rice Narrow leaf1 gene, which encodes a novel protein, affects vein patterning and polar auxin transport. Plant Physiolgy 2008, 147, 1947–1959. [Google Scholar] [CrossRef] [Green Version]
  15. Cho, S.H.; Yoo, S.C.; Zhang, H.T.; Pandeya, D.; Koh, H.J.; Hwang, J.Y.; Kim, G.T.; Paek, N.C. The rice narrow leaf2 and narrow leaf3 loci encode WUSCHEL-related homeobox 3A (OsWOX3A) and function in leaf, spikelet, tiller, and lateral root development. New Phytol. 2013, 198, 1071–1084. [Google Scholar] [CrossRef]
  16. Fujino, K.; Matsuda, Y.; Ozawa, K.; Nishimura, T.; Koshiba, T.; Fraaije, M.W.; Sekiguchi, H. NARROW LEAF 7 controls leaf shape mediated by auxin in rice. Mol. Genet. Genom. 2008, 279, 499–507. [Google Scholar] [CrossRef] [Green Version]
  17. Hu, J.; Zhu, L.; Zeng, D.L.; Gao, Z.Y.; Guo, L.B.; Fang, Y.X.; Zhang, G.H.; Dong, G.J.; Yan, M.X.; Liu, J.; et al. Identification and characterization of NARROW AND ROLLED LEAF 1, a novel gene regulating leaf morphology and plant architecture in rice. Plant Mol. Biol. 2010, 73, 283–292. [Google Scholar] [CrossRef]
  18. Ma, L.; Sang, X.C.; Zhang, T.; Yu, Z.Y.; Li, Y.F.; Zhao, F.M.; Wang, Z.W.; Wang, Y.T.; Yu, P.; Wang, N.; et al. ABNORMAL VASCULAR BUNDLES regulates cell proliferation and procambium cell establishment during aerial organ development in rice. New Phytol. 2017, 213, 275–286. [Google Scholar] [CrossRef] [Green Version]
  19. Garg, R.; Singh, V.K.; Rajkumar, M.S.; Kumar, V.; Jain, M. Global transcriptome and coexpression network analyses reveal cultivar-specific molecular signatures associated with seed development and seed size/weight determination in chickpea. Plant J. 2017, 91, 1088–1107. [Google Scholar] [CrossRef] [Green Version]
  20. Santosh, K.; Santosh, T.; Prasad, S.S.; Archana, P.; Fahamida, A.; Abu, S.M.; Jyothi, B.; Prasad, D.S.; Rudra, B.; Natividad, M. ARice breeding for yield under drought has selected for longer flag leaves and lower stomatal density. J. Exp. Bot. 2021, 72, 4981–4992. [Google Scholar]
  21. Kusterer, B.; Muminovic, J.; Utz, H.F.; Piepho, H.P.; Barth, S.; Heckenberger, M.; Meyer, R.C.; Altmann, T.; Melchinger, A.E. Analysis of a triple testcross design with recombinant inbred lines reveals a significant role of epistasis in heterosis for biomass-related traits in Arabidopsis. Genetics 2007, 175, 2009–2017. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  22. Chen, F.; Hu, Y.; Alessandro, V.; Wu, K.C.; Cai, H.Y.; Qin, Y.; Alison, M.; Lin, Z.G.; Zhang, L.S. The WRKY transcription factor family in model plants and crops. Crit. Rev. Plant Sci. 2018, 36, 311–335. [Google Scholar] [CrossRef]
  23. Gregersen, P.L.; Holm, P.B.; Krupinska, K. Leaf senescence and nutrient remobilisation in barley and wheat. Plant Biol. 2008, 10, 37–49. [Google Scholar] [CrossRef] [PubMed]
  24. Kong, L.G.; Guo, H.H.; Sun, M.Z. Signal transduction during wheat grain development. Planta 2015, 241, 789–801. [Google Scholar] [CrossRef] [PubMed]
  25. Stigter, K.; Plaxton, W. Molecular Mechanisms of Phosphorus Metabolism and Transport during Leaf Senescence. Plants 2015, 4, 773–798. [Google Scholar] [CrossRef] [Green Version]
  26. Wang, F.M.; Rose, T.; Jeong, K.; Kretzschmar, T.; Wissuwa, M. The knowns and unknowns of phosphorus loading into grains, and implications for phosphorus efficiency in cropping systems. J. Exp. Bot. 2016, 67, 122. [Google Scholar] [CrossRef] [Green Version]
  27. Julia, C.; Wissuwa, M.; Kretzschmar, T.; Jeong, K.; Rose, T. Phosphorus uptake, partitioning and redistribution during grain filling in rice. Ann. Bot. 2016, 118, 1151–1162. [Google Scholar] [CrossRef] [Green Version]
  28. Bian, J.M.; Jiang, L.; Liu, L.L.; Xiao, Y.H.; Wang, Z.Q.; Zhao, Z.G.; Zhai, H.Q.; Wan, J.M. Identification of Japonica Chromosome Segments Associated with Heterosis for Yield in Indica × Japonica Rice Hybrids. Crop Sci. 2012, 50, 2328. [Google Scholar] [CrossRef]
  29. Stuber, C.; Edwards, M.; Wendel, J. Molecular Marker-Facilitated Investigations of Quantitative Trait Loci in Maize.II. Factors Influencing Yield and its Component Traits. Crop Sci. 1987, 27, 639–648. [Google Scholar] [CrossRef] [Green Version]
  30. Rajkumar, A.P.; Qvist, P.; Lazarus, R.; Lescai, F.; Ju, J.; Nyegaard, M.; Mors, O.; Børglum, A.D.; Li, Q.; Christensen, J.H. Experimental validation of methods for differential gene expression analysis and sample pooling in RNA-seq. BMC Genom. 2015, 25, 548. [Google Scholar] [CrossRef] [Green Version]
  31. Audic, S.; Claverie, J. The significance of digital gene expression profiles. Genome Res. 1997, 7, 986. [Google Scholar] [CrossRef] [PubMed]
Plants 12 02496 g001 550
Figure 1. DEGs in super hybrid WFYT025 combination. (A) Number of DEGs between the hybrid and its parents in different stages under different environments. (B) Venn diagram of DEGs between the hybrid and its parents under the early rice environment (C) Venn diagram of DEGs between the hybrid and its parents under the middle rice environment. CH, WFY and WFB represent CHT025, WFYT025 and WFB, respectively. E and M represent early and middle rice environments.
Figure 1. DEGs in super hybrid WFYT025 combination. (A) Number of DEGs between the hybrid and its parents in different stages under different environments. (B) Venn diagram of DEGs between the hybrid and its parents under the early rice environment (C) Venn diagram of DEGs between the hybrid and its parents under the middle rice environment. CH, WFY and WFB represent CHT025, WFYT025 and WFB, respectively. E and M represent early and middle rice environments.
Plants 12 02496 g001
Plants 12 02496 g002 550
Figure 2. Gene ontology (GO) enrichment analysis of DEGs (A) GO enrichment analysis of DEGs at 1DPA and at 10 DPA under the environment of early rice in CHT025. (B) GO enrichment analysis of DEGs at 1DPA and at 10 DPA under the environment of early rice in WFB.
Figure 2. Gene ontology (GO) enrichment analysis of DEGs (A) GO enrichment analysis of DEGs at 1DPA and at 10 DPA under the environment of early rice in CHT025. (B) GO enrichment analysis of DEGs at 1DPA and at 10 DPA under the environment of early rice in WFB.
Plants 12 02496 g002
Plants 12 02496 g003 550
Figure 3. KEGG pathway assignments of DEGs. (A) KEGG analysis of DEGs between CHT025 in the environment of early rice. (B) KEGG analysis of DEGs between WFB in the environment of early rice. (A,B) showed the top 20 most represented categories and the number of transcripts predicted to belong to each category.
Figure 3. KEGG pathway assignments of DEGs. (A) KEGG analysis of DEGs between CHT025 in the environment of early rice. (B) KEGG analysis of DEGs between WFB in the environment of early rice. (A,B) showed the top 20 most represented categories and the number of transcripts predicted to belong to each category.
Plants 12 02496 g003
Plants 12 02496 g004 550
Figure 4. Analysis results of 5000 genes in co-expression module. The top half of the graph is divided into a clustering tree of 5000 genes based on their expression profiles into 36 sets of RNA Seq data, with each branch representing one gene.The bottom half of the graph is divided into co-expression modules corresponding to the 5000 genes, and the same co-expression modules are represented by the same color. The analysis resulted in 24 co-expression modules.
Figure 4. Analysis results of 5000 genes in co-expression module. The top half of the graph is divided into a clustering tree of 5000 genes based on their expression profiles into 36 sets of RNA Seq data, with each branch representing one gene.The bottom half of the graph is divided into co-expression modules corresponding to the 5000 genes, and the same co-expression modules are represented by the same color. The analysis resulted in 24 co-expression modules.
Plants 12 02496 g004
Plants 12 02496 g005 550
Figure 5. Go enrichment analysis of midlightblue module about 106 genes in the co-expression module showed statistically significant GO term (FDR < 0.05). The x-axis indicates the level of statistical significance, and the y-axis indicates the significantly enriched GO term.
Figure 5. Go enrichment analysis of midlightblue module about 106 genes in the co-expression module showed statistically significant GO term (FDR < 0.05). The x-axis indicates the level of statistical significance, and the y-axis indicates the significantly enriched GO term.
Plants 12 02496 g005
Plants 12 02496 g006 550
Figure 6. Breakdown of the SDGhps according to the dominance ratio Hp under the environment of early rice. Depending on the principle of Hp = [d]/[a], Hp was classified as either positive or negative. (A) Breakdown of the SDGhps according to the dominance ratio Hp On the 1st day after DAA. (B) Breakdown of the SDGhps according to the dominance ratio Hp On the 10th day after DAA.
Figure 6. Breakdown of the SDGhps according to the dominance ratio Hp under the environment of early rice. Depending on the principle of Hp = [d]/[a], Hp was classified as either positive or negative. (A) Breakdown of the SDGhps according to the dominance ratio Hp On the 1st day after DAA. (B) Breakdown of the SDGhps according to the dominance ratio Hp On the 10th day after DAA.
Plants 12 02496 g006
Plants 12 02496 g007 550
Figure 7. Breakdown of the SDGhps according to the dominance ratio Hp under the environment of middle rice. Depending on the principle of Hp = [d]/[a], Hp was classified as either positive or negative. (A) Breakdown of the SDGhps according to the dominance ratio Hp on the 1st day after DAA. (B) Breakdown of the SDGhps according to the dominance ratio Hp On the 10th day after DAA.
Figure 7. Breakdown of the SDGhps according to the dominance ratio Hp under the environment of middle rice. Depending on the principle of Hp = [d]/[a], Hp was classified as either positive or negative. (A) Breakdown of the SDGhps according to the dominance ratio Hp on the 1st day after DAA. (B) Breakdown of the SDGhps according to the dominance ratio Hp On the 10th day after DAA.
Plants 12 02496 g007
Table
Table 1. Phenotypic analysis of WFYT025 and its parents of early rice and middle rice.
Table 1. Phenotypic analysis of WFYT025 and its parents of early rice and middle rice.
EnvironmentTraitsWFBCHT025WFYT025MPH (%)HPH (%)
       
early rice1D leaf length (cm)28.17 ± 1.48525.2 ± 1.9230.4 ± 2.3213.9 *7.9
1D leaf width (cm)2.02 ± 0.1072.18 ± 0.332.38 ± 0.3013.3 *9.1 *
1D leaf area (cm2)42.6841.254.2629.3 **27.1 **
1D 1000 spikelet weight (g)6.75.687.317.9 **8.9
       
10D leaf length (cm)31.8 ± 3.58530.15 ± 1.8030.18 ± 1.29−2.6−5
10D leaf width (cm)2.07 ± 0.0802.33 ± 0.222.35 ± 0.286.80.8
10D leaf area (cm2)49.3752.6953.194.20.1
10D TGW (g)24.8614.7720.422.81−17.8 **
       
middle rice1D leaf length (cm)2738.544.134.7 **14.5
1D leaf width (cm)1.9322.16.95
1D leaf area (cm2)39.0857.7569.4643.4 **20.2 **
1D 1000 spikelet weight (g)5.034.254.844.3−3.7
       
10D leaf length (cm)31.3740.7843.6821.1 **7.1
10D leaf width (cm)2.112.222.285.32.7
10D leaf area (cm2)49.6467.974.6927.1 **10
10D TGW (g)29.0614.6421.1−3.4−27.3 **
** Significant difference with p < 0.01, * Significant difference with p < 0.05.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Cheng, Q.; Huang, S.; Lin, L.; Zhong, Q.; Huang, T.; He, H.; Bian, J. Genetic Analysis for the Flag Leaf Heterosis of a Super-Hybrid Rice WFYT025 Combination Using RNA-Seq. Plants 2023, 12, 2496. https://doi.org/10.3390/plants12132496

AMA Style

Cheng Q, Huang S, Lin L, Zhong Q, Huang T, He H, Bian J. Genetic Analysis for the Flag Leaf Heterosis of a Super-Hybrid Rice WFYT025 Combination Using RNA-Seq. Plants. 2023; 12(13):2496. https://doi.org/10.3390/plants12132496

Chicago/Turabian Style

Cheng, Qin, Shiying Huang, Lan Lin, Qi Zhong, Tao Huang, Haohua He, and Jianmin Bian. 2023. "Genetic Analysis for the Flag Leaf Heterosis of a Super-Hybrid Rice WFYT025 Combination Using RNA-Seq" Plants 12, no. 13: 2496. https://doi.org/10.3390/plants12132496

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop