Non-coding RNA expression analysis revealed the molecular mechanism of flag leaf heterosis in inter-subspecific hybrid rice

Heterosis has been used widespread in agriculture, but its molecular mechanism is inadequately understood. Plants have a large number of non-coding RNAs (ncRNAs), among them, functional ncRNAs that have been studied widely containing long non-coding RNA (lncRNA) and circular RNA (circRNA) that play a role in varied biological processes, as well as microRNA (miRNA), which can not only regulate the post-transcriptional expression of target genes, but also target lncRNA and circRNA then participate the competing endogenous RNA (ceRNA) regulatory network. However, the influence of these three ncRNAs and their regulatory relationships on heterosis is unknown in rice. In this study, the expression profile of ncRNAs and the ncRNA regulatory network related to heterosis were comprehensively analyzed in inter-subspecific hybrid rice. A total of 867 miRNAs, 3,278 lncRNAs and 2,521 circRNAs were identified in the hybrid and its parents. Analysis of the global profiles of these three types of ncRNAs indicated that significant differences existed in the distribution and sequence characteristics of the corresponding genes. The numbers of miRNA and lncRNA in hybrid were higher than those in its parents. A total of 784 ncRNAs (169 miRNAs, 573 lncRNAs and 42 circRNAs) showed differentially expressed in the hybrid, and their target/host genes were vital in stress tolerance, growth and development in rice. These discoveries suggested that the expression plasticity of ncRNA has an important role of inter-subspecific hybrid rice heterosis. It is worth mentioning that miRNAs exhibited substantially more variations between hybrid and parents compared with observed variation for lncRNA and circRNA. Non-additive expression ncRNAs and allele-specific expression genes-related ncRNAs in hybrid were provided in this study, and multiple sets of ncRNA regulatory networks closely related to heterosis were obtained. Meanwhile, heterosis-related regulatory networks of ceRNA (lncRNA and circRNA) and miRNA were also demonstrated.

Introduction

Heterosis is a ubiquitous phenomenon in sexually reproducing organisms, which has been widely used in crops, and has achieved great social and economic benefits. To date, three genetic hypotheses for heterosis were proposed, which included dominance, overdominance and epistasis (Larièpe et al., 2012; Shen et al., 2014). Due to the complex molecular mechanisms, the current explanations for heterosis are still controversial (Xiong et al., 2021).

The yield of hybrid rice was 20–30% higher than conventional rice (Qian et al., 2021). The research on heterosis of hybrid rice is crucial for improving rice yield and has been the focus of scientists for many years (Shao et al., 2019; Xiong et al., 2021). Two traditional subspecies, indica (Oryza sativa L. subsp. indica Kato) and japonica (O. sativa L. subsp. japonica Kato), have different physiological traits, genomic and epigenetic differentiation (Liu et al., 2007; Lu et al., 2022). These differences lead to obvious heterosis in inter-subspecific hybrid rice, which has potential production and research value (Zhang, 2022).

With the rapid rise of modern genomics techniques, noncoding RNA or nonprotein-coding RNAs (ncRNA or npcRNA), once called genomic “dark matter,” has been discovered (Morris and Mattick, 2014; Yu et al., 2019). By targeting to corresponding DNA, RNA and proteins, ncRNAs participate in a variety of biological processes and molecular mechanisms in plants (Ohtani, 2017; Chand Jha et al., 2021). In rice, many experimental reports have indicated the importance of ncRNAs in growth, development, environmental response (Yuan et al., 2018; Luan et al., 2019; Wang et al., 2020; Zhai et al., 2020; Qiao et al., 2021). Regrettably, the molecular functions of ncRNAs in plant heterosis formation remain poorly understood.

According to the length and function, ncRNAs in plants can be into lncRNAs, circRNAs, miRNAs, ribosomal RNAs (rRNAs), transfer RNAs (tRNAs) and others. MiRNA is a kind of single-stranded endogenous ncRNA abundant in animal and plant kingdom, about 20–24 nt in length, which regulates gene expression by pairing with cognate mRNA bases to silenced mRNA (Voinnet, 2009). LncRNA is an ncRNA with a length of more than 200 nt, located in the nucleus and cytoplasm, which can be categorized into long intergenic ncRNA (lincRNA), intron ncRNA (incRNA) and natural antisense transcript (NAT). It is involved in molecular functions such as chromatin level modification, alternative splicing of mRNA and the regulation of post-translation (Bardou et al., 2014; Laurent et al., 2015; Budak et al., 2020). CircRNA is a kind of closed-loop single-stranded ncRNA produced by reverse splicing, which is associated with a variety of biological processes in eukaryotes (Guo et al., 2014; Wang et al., 2019; Zhang et al., 2020). Although more and more circRNA profiles have been analyzed in plants, the role of circRNAs in plants remains limited compared to the numerous advances made in animals (Lu et al., 2015; Dou et al., 2017; Pan et al., 2018; Han et al., 2020).

In addition to the regulation of target mRNA by miRNAs mentioned above, some lncRNAs and circRNAs can competitively bind miRNAs then influence the expression of target mRNAs, which were called endogenous target mimics (eTMs) or competing endogenous RNAs (ceRNAs). Osa-miR2118 triggers the cleavage of lncRNA, Photoperiod sensitive male sterility 1 (PMS1T), and this interaction of miRNA-lncRNA can result in male sterility of photoperiod sensitive rice (Fan et al., 2016). Two lncRNAs, MSTRG.66289.1 and MSTRG.52515.5, have been experimentally verified as ceRNAs of osa-miR156 and osa-miR396 in rice, respectively, which involved in the fertility conversion (Wang et al., 2021c). A total of 4,517 pairs of circRNA/lncRNA–miRNA–mRNA interaction networks were studied in MH63, and these networks might be involved in growth, development, metabolism and so on (Zhou et al., 2021). Despite these prior miRNA–lncRNA or miRNA–circRNA network investigations, we still have an incomplete understand of ceRNA in inter-subspecific hybrid rice.

In this study, we provided a deep exploration of ncRNA expression profiles and inheritance patterns among an inter-subspecific hybrid rice, covering expression analysis of three types of ncRNAs, identification of heterosis-related ncRNAs, and regulation analysis of heterosis-related miRNAs on two other ncRNAs. The additive expression of heterosis-related ncRNAs was analyzed, and ceRNA networks of these heterosis-associated ncRNAs were also demonstrated. This research provided a more detailed view of ncRNAs and their contribution to heterosis.

Materials and methods

Material planting and collection

The inter-subspecific hybrid rice combinations used in this study were proved to have flag leaf heterosis in the hybrid ZY19 (Wang and Wang, 2022). ZY19 and its parents (Z04A, japonica maternal line; ZHF1015, indica paternal line) were planted at test field of Wuhan University, Wuhan, China. Fertilization and pest control were consistent with conventional field management. At heading stage, flag leaves of hybrid and its parents (three biological replicates of each variety) were collected, frozen in liquid nitrogen for 5–6 h and stored at −80°C for RNA extraction.

RNA extraction and quality detection

RNAs from the 9 sample of flag leaves were extracted by RN40 (Aidlab Biotechnologies, Beijing, China). NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, United States) and Agilent Bioanalyzer 2,100 (Agilent Technologies, Santa Clara, CA, United States) detected RNA concentration, purity, and integrity. High-quality RNAs (concentration ≥ 313.9 ng/μl, RNA integrity number > 7.2) were used to prepare libraries needed for further sequencing.

Preparation of rRNA-depleted strand-specific RNA (ssRNA) library and sequencing

The rRNA-depleted ssRNA libraries were used to identify mRNAs, lncRNAs, and circRNAs. The extracted high-quality RNAs were used as materials, and Ribo-off^® rRNA depletion kit (Vazyme, Nanjing, China) was used to remove rRNA. The first-strand cDNA was synthesized after fragmenting the rRNA-depleted RNA. Then, the second-strand cDNA was synthesized and purified. After purification, base A and adaptor were added to the 3′ end. By PCR amplification and quality detection, nine libraries were constructed for sequencing, which was performed on Illumina NovaSeq 6,000 Platform (Illumina, San Diego, CA, United States) to generate pair-ended 150 bp reads.

Transcriptome assembly and lncRNA, circRNA identification

Prior to assembly, low-quality reads were removed from the raw reads of fastq format. After sequencing quality control, a total of 64.14 Gb (Q30 > 94.15%) of clean data were obtained. Clean reads were mapped onto Nipponbare (MSU Rice Genome Annotation Project Release 7, MSU-V7.0)¹ by Hisat2. Finally, String Tie was used to compare reads on pairs for transcription splicing and quantification.

Transcripts annotated with Gffcompare were obtained, and unknown transcripts with length greater than 200 nt and exon number greater than two were used to screen for putative lncRNAs. From these transcripts, Coding Potential Calculator (CPC), Coding-Non-Coding Index (CNCI), Protein family (Pfam), and Coding Potential Assessment Tool (CPAT) were further used to obtain RNA with no coding ability. As well as the different types of lncRNAs including lincRNA, antisense lncRNA, sense lncRNA, and intronic lncRNA were distinguished. Based on the interaction mechanism between lncRNAs and their target genes, there were two methods to predict target genes of lncRNAs. On the one hand, since lncRNA can regulate the expression of their adjacent mRNA, the mRNAs within 100 kb of corresponding lncRNAs were designated as cis-target genes. On the other hand, lncRNA interacts with target mRNA through base complementation and predicted by LncTar.

The data obtained by ssRNA sequencing were used to predict circRNA by Findcirc software. Since the circRNA splicing site could not be directly mapped onto the reference genome, 20 bp at each end of unmapped reads was taken as anchor, which will be compared to the genome as an independent unit. After finding the binding position of circRNAs, if the sequences on both sides were GT/AG, they were judged to be circRNAs. Candidate circRNA screening pipeline was developed as follows. (1) GU/AG appeared on both sides of the splicing site. (2) Specific breakpoints. (3) Only two mismatches. (4) Maximum length of breakpoint was 2 nt. (5) The corresponding junction was supported by two reads at least. (6) Mapping onto the correct short sequence scored more than 35 points higher than mapping onto other positions. Finally, candidate circRNAs were identified.

Small RNA library construction and sequencing

Using total RNAs extracted as input materials, nine sequencing libraries were generated using the VAHTS Small RNA Library Prep Kit for Illumina (Vazyme, Nanjing, China). The 3′ and 5′ ends of small RNA were connected to universal adapters, respectively, and then were amplified by reverse transcription and PCR. The target fragments were screened by gel separation techniques, and small RNA libraries suitable for Illumina platform were obtained. NovaSeq6000 (Illumina, San Diego, CA, United States) was used for small RNA sequencing after library quality inspection.

miRNA identification and target prediction

The small RNA sequencing original data quality control process was as follows: (1) Removal of the joint. (2) Removal of sequences shorter than 18 nt or longer than 30 nt. (3) Removal of sequences with low-quality values for each sample. (4) Remove reads with N (N is an unidentifiable base) greater than 10%. Using Bowtie software, the clean reads obtained were compared with Silva database, GtRNAdb database, Rfam database and Repbase database for short sequences. SnoRNA, rRNA, tRNA and repeat sequences were filtered out to obtain unannotated reads containing miRNA. Then these obtained reads were mapped onto Nipponbare (MSU-V7.0) by Bowtie and obtained mapped reads. Then the sequences were compared with miRBase (V22) to identify the known miRNA. The novel miRNAs were predicted by miRDeep2. TargetFinder predicted miRNA targets based on known and novel miRNA and corresponding species genomic information.

Differential expression and non-additive expression analysis

Differential expression analysis of four types of RNAs in three varieties (three biological replicates per variety) was analyzed by DESeq2. Fold Change (FC) represents the expression ratio between two varieties and p value of the original hypothesis represents the probability of expressing indifference. The identification of differentially expressed genes (DEGs), differentially expressed lncRNAs (DElncRNAs), differentially expressed circRNAs (DEcircRNAs) and differentially expressed miRNAs (DEmiRNAs) was based on the following thresholds: $\left|{log}_2\mathrm{fold\; change}\right|>1$ and p < 0.05. The differentially expressed RNAs between hybrid and parents (including DEmRNA_HP, DEmiRNA_HP, DElncRNA_HP, and DEcircRNA_HP) were speculated to be involved in heterosis. Correspondingly, the differentially expressed RNAs between two parents were represented by DEmRNA_PP, DEmiRNA_PP, DElncRNA_PP, and DEcircRNA_PP.

Additive and non-additive ncRNA expression analysis was performed by comparing the ncRNA expression levels of hybrids with mid-parent value (MPV). The thresholds to identified non-additive ncRNAs were:

$\left|{log}_2\left(\begin{array}{l}\mathrm{ncRNA\; expression\; of} MPV\\ /\mathrm{ncRNA\; expression\; in} F1 \mathrm{hybrid}\end{array}\right)\right|>1 \mathrm{and} p < 0.05.$

Heterosis-related ncRNA–mRNA and ncRNA–ncRNA integration analysis

To obtain heterosis-related ncRNA-mRNA regulatory network, we analyzed the regulatory relationship between ncRNAs and their corresponding targeted allele-specific expression genes (ASEGs) in hybrid. The ASEGs data are derived from previous studies (Wang and Wang, 2022) of the same inter-subspecific hybrid rice in this study. DEmRNA_HPs corresponding to DEmiRNA_HPs and DElncRNA_HPs were analyzed. Meanwhile, DElncRNA_HPs and DEcircRNA_HPs regulated by DEmiRNA_HPs were also analyzed. These networks were visualized using Cytoscape.

Using weighted correlation network analysis (WGCNA) in R software, the correlation coefficient between non-additive expressed ncRNAs and DEmRNA_HPs was evaluated. Correlation coefficient with corresponding > 0.8 (positive) or < 0.8 (negative) was considered to be a co-expressed mRNA of non-additive expressed ncRNAs.

ncRNA expression level validation

The results of the ssRNA and miRNA sequencing were validated by quantitative real-time PCR (qRT-PCR). A total of 10 lncRNAs, 5 circRNAs, and 5 miRNAs were randomly selected for verification. Total RNA was extracted from nine flag Leaf samples (including three biological replicates) using TRIzol reagent (Catalog No.15596–026, Invitrogen, Carlsbad, CA, United States). LncRNAs and circRNAs were verified using cDNA obtained by reverse transcription with random primers. According to the sequence of miRNA, cDNA synthesized with reverse transcription primer with stem-loop was used for quantitative verification of miRNA. SYBR (PowerUp^™ SYBR^™ Green Master Mix, Thermo Scientific^™, United States) were used on an ABI Step One Plus Real-Time PCR System. OsActin1 was a housekeeping gene for quantitative validation of lncRNA and circRNA. Forward primers and reverse universal primers for miRNA were designed according to the sequences of miRNA and reverse transcription primers, and U6 was used as housekeeping gene. The results were analyzed by comparative Ct (2^–ΔΔCt).

Results

ncRNA and mRNA characteristics in the inter-subspecific hybrid rice combination

A total of 110 million clean reads were obtained from 9 small RNA libraries (Supplementary Table 1). Through mapping to the miRBase (V22) database and prediction by miRDeep2, 867 miRNAs were identified, containing 345 novel miRNAs and 522 known miRNAs (Table 1). Among these identified miRNAs, 514 miRNAs belonged to 153 gene families. To profile lncRNAs, circRNA and mRNA in each sample, ssRNA libraries were successfully constructed using the same plant material used for miRNA sequencing. After filtering low-quality reads, an average of 47.96 million reads were generated in each sample (Supplementary Table 1), and 88.98–91.90% of the clean reads were mapped uniquely to a single best location in the reference genome. Then a total of 3,278 lncRNAs were identified through basic screening and coding ability screening, including lincRNA 1789 (54.58%), antisense lncRNA 955 (29.13%), sense lncRNA 315 (9.61%), and intronic lncRNA 219 (6.68%; Table 1). Find_circ was used to predict the total circular RNA of three varieties, and 2,521 circRNAs were identified for subsequent analysis (Table 1). The mRNA expression levels were also profiled by lncRNA sequencing and the number of identified mRNAs was shown in Table 1. Furthermore, qRT-PCR was used to verify the relative expression levels of the three types of ncRNAs identified among different varieties, which indicated that the sequencing results were reliable (Supplementary Figure 1).

TABLE 1

Table 1. Number of RNA identified in the inter-subspecific hybrid rice combination.

To investigate the differences in chromosome distribution, sequence and structure of different types of RNAs corresponding genes, the characteristics of the corresponding DNA sequences of the identified RNAs were compared in this study. Analysis of genomic characteristics of corresponding genes of RNAs showed that the proportion of all identified RNA corresponding sequences on chromosome 1 was higher than that on other chromosomes, followed by chromosome 2, and this distribution pattern might be related to chromosome length (Supplementary Figure 2A). Density statistics of all identified mRNA, miRNA, lncRNA, and circRNA corresponding genes were carried out with 10,000 bp sliding window. The results showed that the corresponding genes of the four types of RNAs were distributed in low density near the centromere on each chromosome (Figure 1A). The density distribution of miRNA and circRNA corresponding genes on each chromosome was similar, but there were many differences in the distribution of miRNA and lncRNA genes on each chromosome (Figure 1A). Furthermore, the proportion of genes corresponding to different types of RNA on each chromosome was calculated. mRNAs corresponding genes accounted for the highest proportion on each chromosome, followed by lncRNAs and circRNAs (Supplementary Figure 2B).

FIGURE 1

Figure 1. Genomic characteristics and structure characteristics of all identified RNAs corresponding genes. (A) Density map of mRNA, miRNA, lncRNA and circRNA corresponding genes on different chromosomes. (B) Sequence length of identified ncRNAs and mRNAs. (C) Corresponding DNA sequence ORF length of identified ncRNAs and mRNAs. (D) The number of TE contained in the corresponding DNA sequence of ncRNAs and mRNAs. (E) The number of Exon contained in the corresponding DNA sequence identified ncRNAs and mRNAs.

The transcript length of miRNA fragments was 18–25 nt, and the number of miRNAs with transcript lengths of 21 and 24 nt was the majority (Figure 1B). Compared with mRNAs, the transcript length of lncRNAs and circRNAs was mainly distributed in the range of 200–400 nt, which was much smaller than that of mRNA. This characteristic was consistent with previous studies on lncRNA in young panicles of rice (Wang et al., 2021c). At the same time, referring to the study on exons, ORF length and TE number of genes corresponding to lncRNAs in young panicles (Wang et al., 2021c), these elements of genes corresponding to 4 types of RNAs were analyzed and counted. MiRNAs corresponding DNA sequences were too short to identify ORF. The ORF length of most lncRNAs and circRNA corresponding genes was less than 200 nt, and 87.5% ORF of the lncRNAs corresponding genes was less than 100 nt, while the ORF length of most mRNAs corresponding genes was greater than 100 nt (Figure 1C). By comparison with the TE region of rice, the number of whole TE contained in RNA corresponding genes was very small (Figure 1D). It was speculated that many RNA corresponding DNA sequences overlap with TE region, but not completely contained. Slightly different from the results of previous studies, most lncRNAs and miRNAs corresponding genes contained no exons, except 57.14% circRNAs corresponding genes (Figure 1E). All mRNAs corresponding genes had exons, and most of them had fewer than 10 exons.

The transcription factor binding sites (TFBSs) of all detected ncRNA and mRNA corresponding genes loci were predicted by TSPTFBS (Liu et al., 2021). The results showed that TFBSs in mRNA corresponding genes loci were significantly more abundant than those in ncRNA corresponding genes loci. ANAC017 was the most common binding site of ncRNA corresponding genes loci, except for small RNAs identified in this study (Supplementary Figure 3). MiRNAs corresponding genes loci had significantly fewer binding sites than other ncRNAs.

Referring to the study of Crisp et al. (2020), the genome was divided into six meta-feature categories, including TE, noncoding RNA transcript loci, genic, intergenic, gene-proximal, and no-annotation. According to the ncRNA corresponding DNA sequences, the distribution of ncRNA corresponding genes in each genomic feature was determined (Supplementary Figure 4A). The number of lncRNAs corresponding DNA sequences overlapped with TE regions was higher than other regions. The corresponding DNA sequences of miRNAs (21, 22, and 24 nt) were mainly overlapped in TE, gene proximal and inter-gene regions. Further analysis found that the proportion of circRNAs corresponding DNA sequences located in exons was higher than that of other ncRNAs (Supplementary Figure 4B).

Global ncRNA profiles of the hybrid and its parents

After analyzing the sequence characteristics of all identified RNAs, the differences in ncRNA profiles between the hybrid and its parents were analyzed. Figure 2A shows the expression level of different types of RNA (mRNA, lncRNA, circRNA and miRNA from outside to inside) on different chromosomes of the three varieties, which showed that the expression of different types of RNAs was different.

FIGURE 2

Figure 2. ncRNA profiles of the hybrid and its parents. (A) Location information and expression levels of four types of RNAs in the three varieties. mRNA, lncRNA, circRNA and miRNA from outside to inside. Different colors indicate different varieties. (B) The number of RNA expressed by different varieties. (C) Expression levels of four RNAs in the three varieties. (D) Expression levels of miRNA of different length in three varieties. *** means significant difference with p < 0.001, ** means significant difference with p < 0.01, and * means significant difference with p < 0.05.

Among the 867 identified miRNAs, 814, 817 and 765 miRNAs were expressed in maternal, hybrid and paternal lines, respectively (Figure 2B). The number of miRNAs expressed in all three varieties (co-expressed miRNA) accounted for the highest proportion, reaching 83.62%. The number of miRNAs expressed in both maternal and hybrid lines, but not expressed in paternal line, was the next (6.11%). Specifically expressed miRNAs in maternal (pattern 1), hybrid (pattern 2) and parental (pattern 3) lines were accounted for 3.46, 2.19 and 1.16% of the total identified miRNAs, respectively. Among 7 expression patterns of miRNAs, the number of pattern 4 (miRNAs that showed no expression in the hybrid line, but expressed in the maternal and paternal lines) was the least (Figure 2B).

Among the 3,278 lncRNAs identified in this study, 2,191, 2,501 and 2,229 lncRNAs were expressed in maternal, hybrid and paternal lines, respectively (Figure 2B). The result showed a larger number of lncRNAs in hybrid than that of its parents. The number of co-expressed lncRNAs was the highest (37.61%). Pattern 6 (lncRNAs expressed in both maternal and hybrid lines, but not in paternal lines) and Pattern 5 (lncRNAs expressed in both hybrid and paternal lines) were followed, accounting for 15.86 and 15.16%, respectively. The specific lncRNAs expressed in maternal line (Pattern 1), hybrid line (Pattern 2), and parental line (Pattern 3) accounted for 9.18, 7.66 and 9.64% of the total identified lncRNAs, respectively. Among the 7 lncRNA expression patterns, pattern 4 had the least amount (4.88%; Figure 2B).

In the total of 2,521 circRNAs, 913, 1,059 and 1,106 circRNAs were expressed in maternal, hybrid and paternal lines, respectively (Figure 2B). The proportion of circRNAs specifically expressed in three varieties (Pattern 1, Pattern 2, and Pattern 3) was higher than that of other patterns, 24.7, 27.56 and 30.94%, respectively. The pattern with the lowest proportion was Pattern 4, with just 2.3%.

As for mRNAs, 29,903, 30,450 and 29,091 transcripts were identified in the maternal, hybrid and paternal lines, respectively (Figure 2B). A total of 65.87% mRNA transcripts were expressed in all three varieties (Pattern 7). The proportion of mRNA with other expression patterns ranged from 2.46 to 7.19%, and pattern 4 had the lowest percentage (2.46%).

The proportion of co-expressed (pattern 7) miRNA was greater than that of mRNA and lncRNA. The proportion of co-expressed circRNAs was the lowest. Additionally, the number of expressed genes in the hybrid was much higher than that in the parents. Analyses of RNA expression patterns showed that mRNA, lncRNA and miRNA were relatively conserved, while the number of circRNAs with specific expression was higher in the three varieties.

The circRNA host genes and target genes of lncRNAs and miRNA were predicted (Supplementary Table 2). Based on the above analysis of RNA expression patterns, the function analysis of the corresponding target genes of seven expression patterns of RNA was carried out. Because the proportion of co-expressed (Pattern 7) miRNAs was the largest, the number of GO terms of target genes enriched was the largest (Supplementary Figure 5). Target genes of specifically expressed miRNAs in the maternal line and paternal line were significantly enriched in 7 and 4 GO items, respectively. In the hybrid, target genes of specific expression of miRNAs were mainly enriched in 5 GO items, including two biological processes (DNA-templated transcription, L-phenylalanine catabolic process) and molecular functions such as phenylalanine ammonialyase activity. There was no overlap among the GO items of the three specifically expressed miRNAs target genes. Interestingly, among the seven expression patterns of lncRNA, target genes of the 6 expression patterns except pattern 4 were significantly enriched in GO:0004190 (Supplementary Figure 6). GO:0004190 is a molecular function term related to aspartic-type endopeptidase activity. Specifically expressed lncRNAs in the hybrid were significantly enriched in GO:0046914 and GO:0016779, which are molecular functional items related to transition metal ion binding and nucleotidyltransferase activity, respectively. In three rice varieties, target genes of specifically expressed circRNA were enriched in the same GO items, such as GO:0009941, GO:0009570 and GO:0009535, which were related to the cell components of chloroplasts (Supplementary Figure 7).

To explore whether there were significant differences in the global expression of RNAs in the three varieties, the expression levels of four types of RNAs were analyzed. The results showed that there were significant differences in mRNA expression levels among different varieties (Figure 2C). The global expression level of miRNA was the highest among all ncRNAs, and it was at the level of high parental expression in the hybrid. The expression level of lncRNA in the hybrid was significantly higher than that in the maternal line, and there was no significant difference in the expression level between the hybrid and paternal lines. The global expression level of circRNA was the lowest among all ncRNAs, and there was no significant difference among the three varieties. MiRNAs of different lengths have various functions, so miRNAs were distinguished by length (Figure 2D). The results showed that most miRNAs except 18 to 20 nt miRNAs were at the mid-parent expression level. Among them, for only 24 nt miRNA, the global expression level in the hybrid was significantly different from that of the maternal line, and the global expression level of the hybrid was close to that of the paternal line.

Differential expression ncRNAs (DEncRNAs) between hybrid and parents

To identify heterosis-related ncRNAs in this hybrid rice combination, DEncRNAs were analyzed. The results indicated that the total number of DEncRNAs between the parents and hybrid was less than the number of DEncRNAs between the two parents for all types of ncRNA (Figure 3A). There were 179 DEmiRNA_HPs identified and 237 corresponding DEmiRNA_PPs. In DElncRNA identification, 573 DElncRNA_HPs and 647 DElncRNA_PPs were obtained. The number of DEcircRNA was the lowest, including 42 DEcircRNA_HPs and 58 DEcircRNA_PPs. Among the DEGs, the number of DEGs between the paternal line and hybrid (1305) was higher than that between the maternal line and hybrid (881; Supplementary Figure 8A). Except for 24 nt DEmiRNAs, the number of DEncRNAs between the hybrid and paternal lines was less than the number of DEncRNAs between the hybrid and maternal lines (Figure 3A). For example, there were 86 DEmiRNAs between the paternal line and hybrid and 115 DEmiRNAs between the maternal line and hybrid. In DElncRNA_HPs, the number of DElncRNA_HPs between the paternal line and hybrid was 311, and that between the maternal line and hybrid was 332.

FIGURE 3

Figure 3. Differential expression ncRNAs between the three varieties. (A) Venn diagram of ncRNA of different difference groups. (B) Proportions of ncRNA with different levels of difference.

In all ncRNAs identified, the proportion of DEmiRNA was higher (Figure 3B; Supplementary Figure 8B). This proportion was higher than that of DElncRNAs and DEcircRNAs. Moreover, compared with the fold change of DEncRNAs between the two parents, the fold change of DEncRNAs between the hybrid and parents was lower, which was especially evident in downregulated ncRNAs. After further studying the genetic pattern of DEncRNA_HPs, it was found that the RNA expression level of the hybrid was influenced by the expression level of its parents (Supplementary Figure 8C). For example, most lncRNAs with higher expression levels in the hybrid than in the maternal line overlapped with DElnc_PPRNAs, which showed high expression levels in the paternal line. This indicated that the high expression of these lncRNAs in the hybrid benefited from the high expression in the paternal line (Supplementary Figure 8C). Similarly, the low expression of ncRNA in the hybrid was also affected by the low expression of parental ncRNA. These results might be related to the formation of mid-parental expression and non-additive expression in the hybrid.

Analysis of additive and non-additive expression of ncRNAs in hybrid

Previous non-additive mRNA analysis showed that the number of additive genes in hybrid was more than that of non-additive, and there was a phenomenon of transgressive up/downregulated expression (Wang and Wang, 2022). Therefore, additive and non-additive analyses were performed on the identified ncRNAs in this study. Additive expression analysis showed that the expression level of most ncRNAs in hybrid was close to the mid-parent value except circRNA (Figure 4A). Assessed via Wald tests, ncRNAs in hybrid were divided into 12 expression patterns (Figure 4B). Except for circRNAs, the number of additively expressed ncRNAs was greater than that of non-additively expressed ncRNAs. Equivalent expression, that is, no change in ncRNA expression, accounted for the majority. Meanwhile, the above studies on the expression patterns and differential expression patterns of ncRNAs indicated that most ncRNA were conserved in hybrid and its parents. Therefore, it is necessary to further analyze the non-additive expression of ncRNAs in the hybrid.

FIGURE 4

Figure 4. Additive and non-additive ncRNAs in hybrids. (A) Ratio distribution of hybrid and MPV in different ncRNAs. (B) Additive and non-additive patterns of ncRNAs. The 12 possible differential expression states in the F1 hybrid relative to its diploid parents. Roman numerals indicate the same categorization as used in Yoo et al. (2013).

The differentially expressed genes between hybrid and parents (DEG_HPs) have been found to be involved in the molecular mechanism of heterosis formation (Wei et al., 2009), so the correlation analysis between non-additive expressed ncRNAs and DEG_HPs might be helpful to reveal the effect of ncRNAs on heterosis. WGCNA was used to analyze the regulation of non-additive expressed ncRNAs on DEG_HPs. The results showed that DEG_HPs could be divided into four expression patterns (Figure 5A). Among them, MEgrey represents the genes that were not clustered, and the number was only 24 (1.3% of the total DEG_HPs). As for the 6 non-additive circRNAs (Figure 5B), 2 circRNAs (Chr1:27418212|27,474,925 and Chr6:19549861|19,564,228) were significantly positively correlated with the DEG_HPs expression of MEblue, and Chr4:8301385|8,301,731 was significantly positively correlated with the DEG_HPs expression of MEturquoise module (Figure 5B). Through correlation analysis, there was no significant correlation between the expression of 6 transgressive regulated lncRNA and DEG_HPs in hybrid. Among the 10 lncRNAs in pattern ELD-P, 6 lncRNAs (MSTRG.12597.1, MSTRG.5975.1, MSTRG.29875.1, MSTRG.3452.2, MSTRG.24799.2 and MSTRG.29512.1) were significantly positively correlated with the expression of DEG_HPs in MEturquoise module (Figure 5C). Meanwhile, MSTRG.12597.1 was significantly positively correlated with DEG_HPs expression of MEturquoise module and negatively correlated with 148 MEbrown genes. Most of the lncRNAs in pattern ELD-M (13/14) were significantly positively correlated with the expression of 693 DEG_HPs in MEblue module (Figure 5D).

FIGURE 5

Figure 5. The regulation of non-additive expressed ncRNAs on DEGHPs. (A) Four expression patterns of DEGHPs obtained by WGCNA analysis. (B) Correlation analysis of non-additive expressed circRNAs and DEGHPs. (C) Correlation analysis of ELD-P expressed lncRNAs and DEGHPs. (D) Correlation analysis of ELD-M expressed lncRNAs and DEGHPs.

Heterosis-related lncRNA–mRNA and miRNA–mRNA networks

There are multilevel regulatory relationships between ncRNAs and mRNAs. MiRNAs directly regulate target mRNAs by binding them, and lncRNAs directly regulate mRNAs by cis and trans regulation. Differential expression of RNA between hybrid and parents was considered to be related to heterosis formation. To clarify the regulatory network of heterosis-related RNAs, the network of target genes of DEmiRNA_HPs and DElncRNA_HPs with differential expression between hybrid and its parents was analyzed (Supplementary Figure 9).

According to the expression patterns of DEncRNA_HPs and DEG_HPs in the hybrid, all DERNA_HPs were divided into eight clusters. A total of 413 DEG_HPs had cis targeting relationships with 357 DElncRNA_HPs (Supplementary Figure 9A). The complex DElncRNA_HPs-DEG_HPs cis regulatory networks include one DEG_HP regulated by multiple DElncRNA_HPs, multiple DEG_HPs regulated by one DElncRNA_HPs, and multiple DElncRNA_HPs regulating multiple DEG_HPs. LOC_Os01g27590, which was predicted to be a transposon protein, showed cis relationship with 7 lncRNAs (Figure 6A). Among them, the expression pattern of 3 lncRNAs in the hybrid was the same as that of LOC_Os01g27590, indicating that in hybrid the differential expression of LOC_Os01g27590 might be influenced by the expression of these three lncRNAs. MSTRG.5953.1 could cis target three DEG_HPs simultaneously, and the differential expression level of these three DEGs in hybrids was the same as that of MSTRG.5953.1, which showed lower than that of the maternal line and higher than that of the paternal line (Figure 6A). LOC_Os06g30179, encoding cytochrome P450, was cis-targeted by MSTRG.25125.1 (Figure 6A). Both of them showed a pattern of higher expression in hybrid than in maternal line but lower expression than in paternal line. Therefore, it can be speculated that the expression of LOC_Os06g30179 may be related to the targeting of MSTRG.25125.1. At the same time, there were some DElncRNA_HPs-DEG_HPs networks with the opposite expression pattern to target genes in the hybrid. For example, LOC_Os06g29400 was targeted by both MSTRG.30315.1 and MSTRG.30299.1 (Figure 6A). The expression of MSTRG.30299.1 in the hybrid was higher than that in the maternal line, but lower than that in the paternal line. On the contrary, LOC_Os06g29400 expression level in hybrid was lower than that in maternal line and higher than that in paternal line. Similar opposing regulatory networks also deserve further investigation.

FIGURE 6

Figure 6. Regulatory networks of DElncRNAHPs and DEmiRNAHPs on their target DEGHPs. (A) Regulatory network of DElncRNAHPs on their cis target DEGHPs. (B) Regulatory network of DElncRNAHPs on their trans target DEGHPs. (C) Regulatory network of DEmiRNAHPs on their target DEGHPs. The shapes and colors of the different types of RNA were consistent with those illustrated in Supplementary Figure S9.

A total of 65 DElncRNA_HPs and 34 DEG_HPs had complex trans-targeting relationships (Supplementary Figure 9B). LOC_Os07g23030 was predicted to be a retrotransposon protein and was trans-regulated by 28 DElncRNA_HPs, showed lower expression in the hybrid than in the maternal line, but higher than that in the paternal line (Figure 6B). MSTRG.30299.1 and MSTRG.18052.1 had the opposite expression pattern to LOC_Os07g23030 in hybrid, indicating that this gene might be indirectly involved in the heterosis-related lncRNA–mRNA network. LOC_Os11g12050, a disease resistance-related gene, was trans targeted by MSTRG.16272.3. A similar opposite expression pattern existed in the LOC_Os11g12050—MSTRG.16272.3 network (Figure 6B).

As shown in Supplementary Figure 9C, there were a total of 69 DEmiRNA_HPs corresponding to 91 DEG_HPs, and the relationship was complicated. A complex interaction network of multiple miRNAs targeting multiple identical mRNAs simultaneously emerged. In rice, osa-miR319b can negatively regulate tiller bud elongation, tiller number and yield per plant, which plays a critical role in leaf morphogenesis and organ boundary (Wang et al., 2021b). Osa-miR319b, as a DEmiRNA_HP in the hybrid, was involved in 3 DEG_HPs regulatory networks (Figure 6C). Among them, 2 DEG_HPs (LOC_Os04g24300 and LOC_Os12g12120) were potential disease resistance protein genes. LOC_Os03g57050 and osa-miR319b showed opposite expression patterns in hybrid. The expression level of LOC_Os03g57050 in the hybrid was higher than that in the paternal line, and the expression level of osa-miR319b in the hybrid was lower than that in the paternal. These results suggested that osa-miR319b might be involved in heterosis formation by miRNA-mRNA networks.

Functional annotation analysis of heterosis-related lncRNA–mRNA and miRNA–mRNA network genes

To further investigate the potential role of lncRNA-mRNA and miRNA-mRNA networks in heterosis, the related target genes were annotated. Cis target DEG_HPs of DElncRNA_HPs were enriched to 9 and 11 terms in molecular function and biological process category, respectively (Supplementary Figure 10A). In molecular function, most cis target DEG_HPs were labeled as ADP-binding and enzyme activity-related terms. Three cis-regulated DEG_HPs were significantly enriched (p < 0.01) in terms related to sugar transmembrane transporter activity, and these three genes (LOC_Os05g12320, LOC_Os01g65880 and LOC_Os01g42110) all belonged to MtN3 family. The result indicated that these three genes and their corresponding DElncRNA_HPs were involved in heterosis formation, and might have an effect on monosaccharide or polysaccharide transport in leaf of hybrid. In biological process, stress response terms such as defense response, immune response and response to bacterium were enriched. Cis target DEG_HPs were mainly enriched in KEGG pathways such as cyanoamino acid metabolism, nucleotide or base excision repair (Supplementary Figure 10B).

Different from cis target DEG_HPs, trans target DEG_HPs of DElncRNA_HPs were enriched into fewer molecular function category than those in biological process (Supplementary Figure 11A). There were 11 and 9 GO terms in molecular function and biological process, respectively. Among them, 8 DEG_HPs significantly enriched (p < 0.01) in defense response and response to UV were worthy of attention. KEGG analysis showed that only a few DEG_HPs were enriched, such as other glycan degradation and amino acid metabolism pathways (Supplementary Figure 11B).

Target DEG_HPs of DEmiRNA_HPs were enriched to 10 molecular function, 8 biological process and 2 cellular component categories (Supplementary Figure 12A). In molecular function, most DEG_HPs were labeled to DNA binding and enzyme activity-related terms. Target DEG_HPs were enriched in biosynthetic metabolism of KEGG analysis showed that these heterosis-related DEmiRNA_HPs-DEG_HPs networks were mainly enriched in metabolic pathways (Supplementary Figure 12B).

Interaction network analysis of ASEGs with lncRNA and miRNA in the hybrid

ASEG analysis has been shown to be associated with heterosis, but there are few studies on whether these genes are regulated by ncRNAs (Lv et al., 2019). Analysis of ASEGs previously identified in ZY19 (Wang and Wang, 2022) showed that 775 maternal ASEGs were cis-regulated by 1,294 lncRNAs, showing 2,116 interaction pairs (Supplementary Figure 13). These interaction networks included both complex relationship pairs, such as an ASEG (LOC_Os07g28040, presumably an o-methyltransferase gene) regulated by multiple lncRNAs, and many one-to-one maternal ASEG-lncRNAs relationship pairs. A total of 205 maternal ASEGs were trans-regulated by 207 lncRNAs, presenting 208 pairs of interaction (Supplementary Figure 14). Compared with cis ASEG-lncRNA network, the number of trans ASEG-lncRNAs network was small and the interaction relationship was simple. In ZY19, 512 paternal ASEGs were cis-regulated by 979 lncRNAs, showing 1,382 pairs of interaction (Supplementary Figure 15). A total of 124 genes with paternal expression bias were trans-regulated by 146 lncRNAs, showing 167 pairs of interaction (Supplementary Figure 16). The interaction network between miRNA and ASEG was more complex, and multiple miRNAs participate in the expression of multiple ASEG simultaneously. A total of 379 maternal ASEGs and 386 miRNAs presented 986 pairs of relationship, and 257 paternal ASEGs and 341 miRNAs presented 684 pairs of relationship (Supplementary Figure S17, S18).

circRNA–host gene network analysis and its role in heterosis

Exploring the potential relationship between circRNAs and their host genes is necessary for further study of circRNAs. In order to understand whether circRNA was related to the expression level of host gene, correlation analysis was conducted on the expression level of circRNA (SRPBM) and its parent gene (FPKM) in the three varieties (Figure 7A). The results showed that there was a weak positive correlation between the expression level of circRNA and its host gene, which was consistent with previous research (Zhou et al., 2021).

FIGURE 7

Figure 7. Regulatory networks of circRNAs and their host genes. (A) Correlation analysis of circRNAs and their host genes. (B) Regulatory network of DEcircRNAHPs and their host DEGHPs. The host DEGHPs and DEcircRNAHPs were represented by circle and square. The colors of the different types of RNA were consistent with those illustrated in Supplementary Figure S9. (C) circRNA-host ASEG-M network. The host ASEG-M and circRNA were represented by circle and square. (D) circRNA-host ASEG-P network. The host ASEG-P and circRNA were represented by circle and square.

To clarify potential circRNA–host gene network involved in the formation of heterosis, DEcircRNA_HPs and its corresponding host DEG_HPs were analyzed (Figure 7B). Compared with miRNA and lncRNA, the networks of DEcircRNA_HPs-host DEG_HPs were simple. In 3 groups of DEcircRNA_HPs-host DEG_HPs interaction networks, the expression patterns of DEcircRNA_HPs and its host DEG_HPs in hybrid were the same, which showed that the three circRNA (Chr1:21424262|21,426,198, Chr5:20738743|20,740,923, Chr9:10618503|10,619,042) and their host genes might be involved in inter-subspecific hybrid heterosis. This result might be important for exploring the function and origin of circRNA in heterosis.

As mentioned above, the research on ASEG is of great help in understanding the molecular mechanism of heterosis. Therefore, the regulatory network of host ASEGs and circRNAs in hybrid was worth studying. Among the maternal ASEGs of hybrid, there were 55 host genes corresponding to 67 circRNAs, forming 67 pairs of circRNA–host ASEG network (Figure 7C). Most of them were one-to-one simple networks, but there were also complex networks, such as the network with LOC_Os12g08810 as the core host ASEG, which is also called OsGGP, encodes galactose phosphorylase and is crucial for plant stress resistance and development. It is the host gene of four circRNAs (Chr12:4533082|4,533,439, Chr12:4533962|4,534,779, Chr12:4534992|4,536,869, Chr12:4535025|4,536,869). A total of 33 paternal host ASEGs were identified, corresponding to 43 circRNAs, forming 43 pairs of circRNA–host ASEG network (Figure 7D). In this network, the most complex one consisted of 5 circRNAs and 1 host ASEG-P (RHMBD6).

Heterosis-related miRNA–RNA network and heterosis-associated ceRNA network construction

Some lncRNAs and circRNAs can act as ceRNAs and competitively bind miRNAs, thus affecting the regulation of miRNAs on target mRNAs. Through TargetFinder analysis, 4,693 pairs of miRNA–lncRNA interaction and 6,258 pairs of miRNA–circRNA interaction were obtained among the identified ncRNAs, including 592 miRNAs, 1,546 lncRNAs and 859 circRNAs (Supplementary Table 3). Heterosis-related ceRNA regulatory networks were established by combining miRNA target genes with DEncRNA_HPs in the hybrid.

In the construction of DElncRNA_HPs-DEmiRNA_HPs-DEG_HPs axes, we obtained complex ceRNA networks, and most of them were centered on miRNA (Figure 8A). There were 52 miRNAs, 39 lncRNAs and 77 DEG_HPs, forming 6 groups of ceRNA interaction networks. In other words, 39 lncRNAs related to heterosis were considered as potential ceRNAs to influence the regulation of miRNAs on their target genes. Based on targeting relationships among different networks, complex ceRNA networks were formed (Figure 8A). As shown in Figure 8A, 15 DElncRNA_HPs simultaneously acted as ceRNAs of osa-miR2103, affecting its regulation of LOC_Os10g27480. LOC_Os10g27480 is located on chromosome 10 of rice and its current function has not been studied. These 17 RNAs were differentially expressed in the hybrid, and also identified as potential ceRNAs. Further study of this network may contribute to a better understanding the contribution of ceRNA to heterosis formation. At the same time, there were also some complex networks of multiple DElncRNA_HPs acting as ceRNAs of multiple DEmiRNA_HPs (including some novel miRNAs identified in this study, such as novel_miR_207, novel_miR_181 and so on). For example, the miR2118 family to which osa-miR2118b belongs is a relatively conserved 22 nt miRNA superfamily in plants, which plays important role in plant growth and development, disease resistance and resistance (Zhang et al., 2022b). Osa-miR2118b was differentially expressed in hybrid and its parents, and two lncRNAs (MSTRG.8205.8 and MSTRG.18758.1) were acted as ceRNAs of osa-miR2118b, affecting the regulation of osa-miR2118b on LOC_Os08g42700 (Figure 8A). LOC_Os08g42700 is presumed to be a resistance protein gene. This is also a set of ceRNA network with potential heterosis that deserve further study.

FIGURE 8

Figure 8. Regulatory networks of DEmiRNAHPs and ceRNAs. (A) Regulatory network of DElncRNAHPs as ceRNA. (B) Regulatory network of DEcircRNAHPs as ceRNA. The target DEGHPs, DElncRNAHPs, DEmiRNAHPs and DEcircRNAHPs were represented by circle, triangle, diamond and square, respectively. The colors of the different types of RNA were consistent with those illustrated in Supplementary Figure S9.

There were 26 miRNAs, 16 circRNAs and 37 DEG_HPs in the DEcircRNA_HPs-DEmiRNA_HPs-DEG_HPs interaction network, forming one complex ceRNA interaction network (Figure 8B). The results showed that Chr2:21854485|21,887,097 was the ceRNA of osa-miR2118n and osa-miR2118b, and influence the regulation of these two miRNAs on LOC_Os08g42700 (Figure 8B). Both osa-miR2118n and osa-miR2118b belong to the miR2118 family, which has been mentioned above. At the same time, Chr1:27418212|27,474,925 can competitively bind with osa-miR820a，osa-miR820b and osa-miR820c, thus affecting the regulation of these three miRNAs on target genes. MiR820 is related to methyltransferase synthesis under pest stress, which affects rice growth and development (Tang and Chu, 2017). Most of the target genes corresponding to osa-miR820a，osa-miR820b and osa-miR820c were speculated to be related to retrotransposon protein-coding genes. These results suggest that these ceRNAs are not only involved in the formation of heterosis, but also have important functions in rice.

Discussion

The molecular mechanism of heterosis is still unclear through extensive studies on genome, transcriptome and epigenetics (Larièpe et al., 2012; Shen et al., 2014; Shao et al., 2019; Wang and Wang, 2022; Zhang, 2022). In recent years, several ncRNAs and their corresponding regulatory networks have been proved to play important roles in plant growth and development (Ohtani, 2017; Yuan et al., 2018; Luan et al., 2019; Wang et al., 2020; Zhai et al., 2020; Chand Jha et al., 2021; Qiao et al., 2021). However, there are few studies on the role of ncRNA in plant heterosis so far, which have been conducted in maize and cabbage (Crisp et al., 2020; Hou et al., 2020; Li et al., 2021). As a hybrid rice with high yield potential, indica × japonica hybrid rice has always been an important material for heterosis research (Xiong et al., 2021; Zhang, 2022). This study demonstrated that ncRNA plays an important role in the heterosis of inter-subspecific hybrid rice by comprehensively analyzing the expression profile of ncRNA and the ncRNA regulatory network related to heterosis.

Specific expression of different types of ncRNA in hybrid

In recent years, studies on plant ncRNAs have mainly focused on lncRNAs, circRNAs and miRNAs (Dou et al., 2017; Wang et al., 2018, 2021c; Li et al., 2020; You et al., 2021; Zhou et al., 2021). Due to the differences in the sources and functions of these three types of ncRNAs, they play a role in various plant growth and development processes (Voinnet, 2009; Tang and Chu, 2017; Zhai et al., 2020; Chand Jha et al., 2021). A total of 3,278 lncRNAs, 2,521 circRNAs and 867 miRNAs (including 345 newly predicted miRNAs) were identified in this study. Different types of ncRNAs corresponded to different length of DNA sequences, which might result in different ORF length, TE numbers and exon numbers of lncRNA, circRNA and miRNA corresponding genes. The majority of lncRNA encoded DNA sequences (99.27%) contained 0 to 100 nt ORFs. This result was consistent with the DNA sequence characteristics of corresponding to lncRNA identified in young panicles of rice (Wang et al., 2021c). At the same time, our study founded that a small number of circRNAs (10.24%) corresponding DNA sequences contained ORFs, while miRNA corresponding genes did not contain ORFs. TFBS prediction results showed that the DNA sequences corresponding to circRNAs had more TFBS than the other two ncRNAs (lncRNA and miRNA). ANAC017 is an endoplasmic reticulum binding transcription factor involved in the response of Arabidopsis cell wall to aluminum stress (Tao et al., 2022). Except for miRNA, ANAC017 was the most common predicted binding site for the other two ncRNA (lncRNA and circRNA) corresponding genes. TFBS prediction of miRNA showed only 4 different transcription factors could bind to miRNA corresponding genes. These results will provide clues for the study of lncRNA, circRNA and miRNA in rice.

Previous studies have shown that hybridization can activate or inhibit gene expression in hybrids (Wei et al., 2009; Lv et al., 2019). In this study, we explored the effect of hybridization on ncRNA expression profiles of hybrid and its parents. The results showed that the expression numbers of lncRNA and miRNA in hybrid were higher than those in parents except for circRNAs. And the hybrid-specific expression of ncRNA (pattern 2) may provide valuable clues to the formation of heterosis. Among these ncRNAs, the specific expression of circRNAs was evident. Then, by GO functional annotation of genes targeted by these specific expression ncRNAs, most of them were related to chlorophyll biosynthesis. These results suggested that these specific expression ncRNAs in the hybrid might be related to flag leaf heterosis of inter-subspecific hybrid rice.

In maize, studies on heterosis-related small RNA of different lengths showed that there were different expression patterns of 21, 22, 23 nt small RNA in different maize samples, and overdominance expression of small RNA were rarely occurred in hybrids. In this study, results indicated significant differences between hybrid and parental miRNA profiles. The expression level of 21, 22, 23 nt miRNA in the hybrid showed no difference with the parents, but the expression level of 24 nt miRNA in the hybrid was significantly different with its parents (Figure 2D). This indicated that miRNAs with different functions have different contributions to genetics, suggesting that we need to conduct further analysis of miRNAs with different lengths.

The important role of lncRNA in inter-subspecific hybrid rice heterosis formation

In plant, lncRNAs play key roles in development, reproduction, responses to light, stress responses and many other processes (Wang et al., 2014, 2020; Zhang et al., 2014; Jha et al., 2020). The expression of lncRNAs in rice shows tissue-specific and stage-specific (Zhang et al., 2014). An endogenous lncRNA, TWISTED LEAF (TL), plays a cis-regulatory role on OsMYB60 in leaf morphological development (Liu et al., 2018). For anther-development, LDMAR regulates photoperiod-sensitive male sterility, an essential component of hybrid rice, and the rice epiallele Epi-sp., which encodes an lncRNA, was found to affect panicle architecture (Ding et al., 2012; Luan et al., 2019). The overexpression of lncRNA LAIR increases rice grain yield and upregulates some neighboring gene expression (Wang et al., 2018). Parental expression level dominance lncRNAs were identified in rice backcross introgression lines (Li et al., 2020), but lncRNA expression profile and the important role of lncRNA in heterosis formation of inter-subspecific hybrid rice need to be studied.

In this study, lncRNA expression profile showed variety-specific expression. At the same time, the DElncRNA_HPs-DEG_HPs cis/trans regulatory network related to heterosis was obtained in the inter-subspecific hybrid rice combinations (Supplementary Figures 9A,B). Cis target gene LOC_Os01g27590 corresponding to the 7 lncRNAs merit further research, especially MSTRG.1441.1, MSTRG.1457.6 and MSTRG.1441.9. The differential expression pattern of LOC_Os01g27590 in hybrid relative to its parents was the same as these three lncRNAs above, indicating that the expression of LOC_Os01g27590 in hybrid was related to these three lncRNAs, which may affect the formation of heterosis in hybrid. DElncRNA_HPs-DEG_HPs networks with opposite expression pattern were also worthy of further study. For example, the network of MSTRG.30315.1/MSTRG.30299.1 and their cis target LOC_Os06g29400. The number of trans regulatory networks of lncRNA is small, but they also provide important clues for heterosis-related networks. LOC_Os07g23030 and its corresponding 28 DElncRNA_HPs formed a network that indirectly affected heterosis. In this network, although LOC_Os07g23030 was differentially expressed in hybrid relative to its parents, its corresponding DElncRNA_HPs differential expression pattern was different. In addition to the above-mentioned multiple lncRNA trans-targeting the same gene, MSTRG.15392.2 can simultaneously target two DEG_HPs and the differential expression pattern of these two genes (LOC_Os11g42070 and LOC_Os04g14204) showed the same with MSTRG.15392.2 in hybrid. These DElncRNA_HPs-DEG_HPs networks might play an important role in heterosis formation.

The potential role of circRNAs in the heterosis

Among the circRNAs identified by Dou et al. (2017) in leaves of Arabidopsis, a certain number of circRNA related genes were enriched in the photosynthetic system. Meanwhile, circRNAs are involved in the response of Arabidopsis to heat stress (Pan et al., 2018). In wheat leaves, 88 candidate circRNAs were identified under dehydration stress, revealing a possible relationship between circRNAs and leaf dehydration response (Wang et al., 2017). Studies on the circRNAs expression profile of maize showed that circRNAs were involved in the response of stress pathway, and their expression diversity was affected by genetic background, which indicated that circRNAs may be involved in various regulatory networks in plant (Han et al., 2020). As for rice, it is necessary to study how circRNAs are affected by the genetic differences between indica and japonica and their roles in the formation of inter-subspecific hybrid rice. In the three varieties in this study, circRNAs were expressed at low levels and with high specificity. A total of five heterosis-related circRNAs and their host gene networks were obtained, among which three groups were highly correlated with heterosis. It is speculated that the differential expression of host gene in hybrid affects the corresponding circRNAs. The interaction analysis between ASEG and circRNAs is helpful for further study of the molecular mechanism of heterosis. In this study, 67 pairs of circRNA-host ASEG networks were obtained, which provided clues for further research on the influence of circRNA in the hybrid on allele-specific expression.

Important contribution of miRNAs and their ceRNAs to heterosis formation

In the past few years, several research groups have reported that some miRNAs are involved in regulating rice growth, development and the morphogenesis of shoot architecture (Tang and Chu, 2017; Qiao et al., 2021). OsASI1 is involved in the processing of a miRNA precursor and regulates the abundance of miRNA, which is essential for rice pollen development and flowering (You et al., 2021). Rice miRNAs can coordinate both immune and agronomic traits. MiR168 has been identified to be related to immune response in rice (Wu et al., 2015). A recent study showed that overexpression of miR168 significantly inhibited the expression of AGO1, resulting in increased sensitivity of rice to Magnaporthe oryzae, delayed growth period, and decreased yield-related traits (Wang et al., 2021a). In this study, osa-miR168a-3p and osa-miR168a-5p were highly expressed in hybrid and its parents. Osa-miR535 positively regulates tillering and negatively regulates yield and disease resistance by targeting OsSPL14 (Zhang et al., 2022a). In this study, osa-miR535-5p was highly expressed in the maternal line, but at a low expression level in both the hybrid and the paternal line. It might be that low expression in the paternal line affected the expression level of osa-miR535-5p in the hybrid.

Analysis of miRNA expression profiles in maize and Chinese cabbage suggested the potential contribution of miRNA-mediated regulatory networks to heterosis (Hou et al., 2020; Li et al., 2021). Regrettably, there are few studies on the molecular mechanism of miRNA-mediated regulatory networks and rice heterosis. In the DElncRNA_HPs-DEmiRNA_HPs-DEG_HPs and DEcircRNA_HPs-DEmiRNA_HPs-DEG_HPs networks obtained in this study, 39 and 16 ceRNAs (lncRNA and circRNA) were predicted, respectively. Among them, MSTRG.8205.8, MSTRG.18758.1 and Chr1:27418212|27,474,925 was forecast for ceRNA associated with heterosis. At the same time, osa-miR820c was proved to be involved not only in the resistance response of rice (Tang and Chu, 2017), but also in the ceRNA network of heterosis formation. Members of osa-miR820 and osa-miR2118 families played important roles in these ceRNA networks, and the results would provide clues for subsequent studies on ceRNA networks and heterosis.

In summary, 867 miRNAs, 3,278 lncRNAs and 2,521 circRNAs were identified in the inter-subspecific hybrid rice combination. Some DElncRNA_HPs-DEG_HPs, DEmiRNA_HPs-DEG_HPs and circRNA-host gene interaction networks were obtained by analyzing the differentially expressed RNAs between hybrid and its parents. For example, the network of seven DElncRNA_HPs cis regulated LOC_Os01g27590, and the network of 28 DElncRNA_HPs simultaneously trans regulated LOC_Os07g23030. The DEG_HPs regulatory networks involved in osa-miR319b, which plays an important role in rice leaves, were also analyzed. This study demonstrated multiple ncRNA regulatory networks in which ASEG participated. These might provide reference networks for further study of the effects of ncRNAs (miRNA, lncRNA, and circRNA) on heterosis. At the same time, we obtained two groups of complex DElncRNA_HPs and DEcircRNA_HPs as ceRNA to involved in the regulatory network of miRNA on target genes. These results not only indicated that the expression plasticity of ncRNAs was important for the heterosis of inter-subspecific hybrid rice, but also provides clues for the further analysis of ncRNA-heterosis-related gene networks.

Data availability statement

Our data have been uploaded to NCBI Sequence Read Archive (SRA) database with accession number SRR19994724-41 (PRJNA855983).

Author contributions

MW performed material preparation, data analysis, and manuscript writing. JW designed the project and modified the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the State Key Basic Research and Development Plan of China (2013CB126900).

Acknowledgments

We are grateful to Linyou Wang (Zhejiang Academy of Agricultural Science, China) for providing rice seed samples.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.990656/full#supplementary-material

Footnotes

1. ^http://rice.plantbiology.msu.edu/

References

Bardou, F., Ariel, F., Simpson, C. G., Romero-Barrios, N., Laporte, P., Balzergue, S., et al. (2014). Long noncoding RNA modulates alternative splicing regulators in Arabidopsis. Dev. Cell 30, 166–176. doi: 10.1016/j.devcel.2014.06.017

PubMed Abstract | CrossRef Full Text | Google Scholar

Budak, H., Kaya, S. B., and Cagirici, H. B. (2020). Long non-coding RNA in plants in the era of reference sequences. Front. Plant Sci. 11:276. doi: 10.3389/fpls.2020.00276

PubMed Abstract | CrossRef Full Text | Google Scholar

Chand Jha, U., Nayyar, H., Mantri, N., and Siddique, K. H. M. (2021). Non-coding RNAs in legumes: their emerging roles in regulating biotic/abiotic stress responses and plant growth and development. Cells 10:1674. doi: 10.3390/cells10071674

PubMed Abstract | CrossRef Full Text | Google Scholar

Crisp, P. A., Hammond, R., Zhou, P., Vaillancourt, B., Lipzen, A., Daum, C., et al. (2020). Variation and inheritance of small RNAs in maize inbreeds and F1 hybrids. Plant Physiol. 182, 318–331. doi: 10.1104/pp.19.00817

PubMed Abstract | CrossRef Full Text | Google Scholar

Ding, J., Lu, Q., Ouyang, Y., Mao, H., Zhang, P., Yao, J., et al. (2012). A long noncoding RNA regulates photoperiod-sensitive male sterility, an essential component of hybrid rice. Proc. Natl. Acad. Sci. U. S. A. 109, 2654–2659. doi: 10.1073/pnas.1121374109

PubMed Abstract | CrossRef Full Text | Google Scholar

Dou, Y., Li, S., Yang, W., Liu, K., Du, Q., Ren, G., et al. (2017). Genome-wide discovery of circular RNAs in the leaf and seedling tissues of Arabidopsis thaliana. Curr. Genomics 18, 360–365. doi: 10.2174/1389202918666170307161124

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, Y., Yang, J., Mathioni, S. M., Yu, J., Shen, J., Yang, X., et al. (2016). PMS1T, producing phased small-interfering RNAs, regulates photoperiod-sensitive male sterility in rice. Proc. Natl. Acad. Sci. U. S. A. 113, 15144–15149. doi: 10.1073/pnas.1619159114

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, J. U., Agarwal, V., Guo, H., and Bartel, D. P. (2014). Expanded identification and characterization of mammalian circular RNAs. Genome Biol. 15:409. doi: 10.1186/s13059-014-0409-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, Y., Li, X., Yan, Y., Duan, M. H., and Xu, J. H. (2020). Identification, characterization, and functional prediction of circular RNAs in maize. Mol. Gen. Genomics. 295, 491–503. doi: 10.1007/s00438-019-01638-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Hou, G., Dong, Y., Zhu, F., Zhao, Q., Li, T., Dou, D., et al. (2020). MicroRNA transcriptomic analysis of the sixth leaf of maize (Zea mays L.) revealed a regulatory mechanism of jointing stage heterosis. BMC Plant Biol. 20:541. doi: 10.1186/s12870-020-02751-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Jha, U. C., Nayyar, H., Jha, R., Khurshid, M., Zhou, M., Mantri, N., et al. (2020). Long non-coding RNAs: emerging players regulating plant abiotic stress response and adaptation. BMC Plant Biol. 20:466. doi: 10.1186/s12870-020-02595-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Larièpe, A., Mangin, B., Jasson, S., Combes, V., Dumas, F., Jamin, P., et al. (2012). The genetic basis of heterosis: multi parental quantitative trait loci mapping reveals contrasted levels of apparent overdominance among traits of agronomical interest in maize (Zea mays L.). Genetics 190, 795–811. doi: 10.1534/genetics.111.133447

PubMed Abstract | CrossRef Full Text | Google Scholar

Laurent, G. S., Wahlestedt, C., and Kapranov, P. (2015). The landscape of long noncoding RNA classification. Trends Genet. 31, 239–251. doi: 10.1016/j.tig.2015.03.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, M., Cao, A., Wang, R., Li, Z., Li, S., and Wang, J. (2020). Genome-wide identification and integrated analysis of lncRNAs in rice backcross introgression lines (BC2F12). BMC Plant Biol. 20:300. doi: 10.1186/s12870-020-02508-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, P., Su, T., Zhang, D., Wang, W., Xin, X., Yu, Y., et al. (2021). Genome-wide analysis of changes in miRNA and target gene expression reveals key roles in heterosis for Chinese cabbage biomass. Hortic Res. 8:39. doi: 10.1038/s41438-021-00474-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, X., Li, D., Zhang, D., Yin, D., Zhao, Y., Ji, C., et al. (2018). A novel antisense long noncoding RNA, TWISTED LEAF, maintains LEAF blade flattening by regulating its associated sense R2R3-MYB gene in rice. New Phytol. 218, 774–788. doi: 10.1111/nph.15023

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, X. H., Lu, T. T., Yu, S. L., Li, Y., Huang, Y. C., Huang, T., et al. (2007). A collection of 10,096 indica rice full-length cDNAs reveals highly expressed sequence divergence between Oryza sativa indica and japonica subspecies. Plant Mol. Biol. 65, 403–415. doi: 10.1007/s11103-007-9174-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, L., Zhang, G., He, S., and Hu, X. (2021). TSPTFBS: a docker image for trans-species prediction of transcription factor binding sites in plants. Bioinformatics 37, 260–262. doi: 10.1093/bioinformatics/btaa1100

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, T., Cui, L., Zhou, Y., Zhu, C., Fan, D., Gong, H., et al. (2015). Transcriptome-wide investigation of circular RNAs in rice. RNA 21, 2076–2087. doi: 10.1261/rna.052282.115

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, Y., Xu, Y., and Li, N. (2022). Early domestication history of Asian rice revealed by mutations and genome-wide analysis of gene genealogies. Rice 15:11. doi: 10.1186/s12284-022-00556-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Luan, X., Liu, S., Ke, S., Dai, H., Xie, X. M., Hsieh, T. F., et al. (2019). Epigenetic modification of ESP, encoding a putative long noncoding RNA, affects panicle architecture in rice. Rice 12:20. doi: 10.1186/s12284-019-0282-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Lv, Z., Zhang, W., Wu, Y., Huang, S., Zhou, Y., Zhang, A., et al. (2019). Extensive allele-level remodeling of histone methylation modification in reciprocal F1 hybrids of rice subspecies. Plant J. 97, 571–586. doi: 10.1111/tpj.14143

PubMed Abstract | CrossRef Full Text | Google Scholar

Morris, K. V., and Mattick, J. S. (2014). The rise of regulatory RNA. Nat. Rev. Genet. 15, 423–437. doi: 10.1038/nrg3722

PubMed Abstract | CrossRef Full Text | Google Scholar

Ohtani, M. (2017). Expanding the plant non-coding RNA world. J. Plant Res. 130, 3–5. doi: 10.1007/s10265-016-0896-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Pan, T., Sun, X., Liu, Y., Li, H., Deng, G., Lin, H., et al. (2018). Heat stress alters genome-wide profiles of circular RNAs in Arabidopsis. Plant Mol. Biol. 96, 217–229. doi: 10.1007/s11103-017-0684-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Qian, Q., Zhang, F., and Xin, Y. (2021). Yuan Longping and hybrid rice research. Rice 14:101. doi: 10.1186/s12284-021-00542-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiao, J., Jiang, H., Lin, Y., Shang, L., Wang, M., Li, D., et al. (2021). A novel miR167a-OsARF6-OsAUX3 module regulates grain length and weight in rice. Mol. Plant 14, 1683–1698. doi: 10.1016/j.molp.2021.06.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Shao, L., Xing, F., Xu, C., Zhang, Q., Che, J., Wang, X., et al. (2019). Patterns of genome-wide allele-specific expression in hybrid rice and the implications on the genetic basis of heterosis. Proc. Natl. Acad. Sci. U. S. A. 116, 5653–5658. doi: 10.1073/pnas.1820513116

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, G., Zhan, W., Chen, H., and Xing, Y. (2014). Dominance and epistasis are the main contributors to heterosis for plant height in rice. Plant Sci. 215-216, 11–18. doi: 10.1016/j.plantsci.2013.10.004

CrossRef Full Text | Google Scholar

Tang, J., and Chu, C. (2017). MicroRNAs in crop improvement: fine-tuners for complex traits. Nat. Plants 3:17077. doi: 10.1038/nplants.2017.77

PubMed Abstract | CrossRef Full Text | Google Scholar

Tao, Y., Wan, J. X., Liu, Y. S., Yang, X. Z., Shen, R. F., and Zhu, X. F. (2022). The NAC transcription factor ANAC017 regulates aluminum tolerance by regulating the cell wall-modifying genes. Plant Physiol. 189, 2517–2534. doi: 10.1093/plphys/kiac197

CrossRef Full Text | Google Scholar

Voinnet, O. (2009). Origin, biogenesis, and activity of plant microRNAs. Cells 136, 669–687. doi: 10.1016/j.cell.2009.01.046

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Chung, P. J., Liu, J., Jang, I. C., Kean, M. J., Xu, J., et al. (2014). Genome-wide identification of long noncoding natural antisense transcripts and their responses to light in Arabidopsis. Genome Res. 24, 444–453. doi: 10.1101/gr.165555.113

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L. L., Jin, J. J., Li, L. H., and Qu, S. H. (2020). Long non-coding RNAs responsive to blast fungus infection in rice. Rice 13:77. doi: 10.1186/s12284-020-00437-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Li, Y., Chern, M., Zhu, Y., Zhang, L. L., Lu, J. H., et al. (2021a). Suppression of rice miR168 improves yield, flowering time and immunity. Nat. Plants 7, 129–136. doi: 10.1038/s41477-021-00852-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Luo, X., Sun, F., Hu, J., Zha, X., Su, W., et al. (2018). Overexpressing lncRNA LAIR increases grain yield and regulates neighbouring gene cluster expression in rice. Nat. Commun. 9:3516. doi: 10.1038/s41467-018-05829-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, M., and Wang, J. (2022). Transcriptome and DNA methylome analyses provide insight into the heterosis in flag leaf of inter-subspecific hybrid rice. Plant Mol. Biol. 108, 105–125. doi: 10.1007/s11103-021-01228-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Xiong, Z., Li, Q., Sun, Y., Jin, J., Chen, H., et al. (2019). Circular RNA profiling of the rice photo-thermosensitive genic male sterile line Wuxiang S reveals circRNA involved in the fertility transition. BMC Plant Biol. 19:340. doi: 10.1186/s12870-019-1944-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, R., Yang, X., Guo, S., Wang, Z., Zhang, Z., and Fang, Z. (2021b). MiR319-targeted OsTCP21 and OsGAmyb regulate tillering and grain yield in rice. J. Integr. Plant Biol. 63, 1260–1272. doi: 10.1111/jipb.13097

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Yang, M., Wei, S., Qin, F., Zhao, H., and Suo, B. (2017). Identification of circular RNAs and their targets in leaves of Triticum aestivum L. under dehydration stress. Front. Plant Sci. 7:2024. doi: 10.3389/fpls.2016.02024

CrossRef Full Text | Google Scholar

Wang, Y., Zhang, H., Li, Q., Jin, J., Chen, H., Zou, Y., et al. (2021c). Genome-wide identification of lncRNAs involved in fertility transition in the photo-thermosensitive genic male sterile rice line Wuxiang S. Front. Plant Sci. 11:580050. doi: 10.3389/fpls.2020.580050

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, G., Tao, Y., Liu, G., Chen, C., Luo, R., Xia, H., et al. (2009). A transcriptomic analysis of superhybrid rice LYP9 and its parents. Proc. Natl. Acad. Sci. U. S. A. 106, 7695–7701. doi: 10.1073/pnas.0902340106

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, J., Yang, Z., Wang, Y., Zheng, L., Ye, R., Ji, Y., et al. (2015). Viral-inducible Argonaute18 confers broad-spectrum virus resistance in rice by sequestering a host microRNA. eLife 4:e05733. doi: 10.7554/eLife.05733

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiong, Y., Zhang, C., Zhou, H., Sun, W., Wang, P., Wang, D., et al. (2021). Identification of heterotic loci with desirable allelic interaction to increase yield in rice. Rice 14:97. doi: 10.1186/s12284-021-00539-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoo, M. J., Szadkowski, E., and Wendel, J. F. (2013). Homoeolog expression bias and expression level dominance in allopolyploid cotton. Heredity (Edinb). 110, 171–180. doi: 10.1038/hdy.2012.94

PubMed Abstract | CrossRef Full Text | Google Scholar

You, L. Y., Lin, J., Xu, H. W., Chen, C. X., Chen, J. Y., Zhang, J., et al. (2021). Intragenic heterochromatin-mediated alternative polyadenylation modulates miRNA and pollen development in rice. New Phytol. 232, 835–852. doi: 10.1111/nph.17635

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, Y., Zhang, Y., Chen, X., and Chen, Y. (2019). Plant noncoding RNAs: hidden players in development and stress responses. Annu. Rev. Cell Dev. Biol. 35, 407–431. doi: 10.1146/annurev-cellbio-100818-125218

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, J., Li, J., Yang, Y., Tan, C., Zhu, Y., Hu, L., et al. (2018). Stress-responsive regulation of long non-coding RNA polyadenylation in Oryza sativa. Plant J. 93, 814–827. doi: 10.1111/tpj.13804

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhai, R., Ye, S., Zhu, G., Lu, Y., Ye, J., Yu, F., et al. (2020). Identification and integrated analysis of glyphosate stress-responsive microRNAs, lncRNAs, and mRNAs in rice using genome-wide high-throughput sequencing. BMC Genom. 21:238. doi: 10.1186/s12864-020-6637-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, G. (2022). The next generation of rice: inter-subspecific Indica-japonica hybrid rice. Front. Plant Sci. 13:857896. doi: 10.3389/fpls.2022.857896

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, L. L., Huang, Y. Y., Zheng, Y. P., Liu, X. X., Zhou, S. X., Yang, X. M., et al. (2022a). Osa-miR535 targets SQUAMOSA promoter binding protein-like 4 to regulate blast disease resistance in rice. Plant J. 110, 166–178. doi: 10.1111/tpj.15663

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, P., Li, S., and Chen, M. (2020). Characterization and function of circular RNAs in plants. Front. Mol. Biosci. 7:91. doi: 10.3389/fmolb.2020.00091

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y. C., Liao, J. Y., Li, Z. Y., Yu, Y., Zhang, J. P., Li, Q. F., et al. (2014). Genome-wide screening and functional analysis identify a large number of long noncoding RNAs involved in the sexual reproduction of rice. Genome Biol. 15:512. doi: 10.1186/s13059-014-0512-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Waseem, M., Zeng, Z., Xu, J., Chen, C., Liu, Y., et al. (2022b). MicroRNA482/2118, a miRNA superfamily essential for both disease resistance and plant development. New Phytol. 233, 2047–2057. doi: 10.1111/nph.17853

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, R., Sanz-Jimenez, P., Zhu, X. T., Feng, J. W., Shao, L., Song, J. M., et al. (2021). Analysis of rice transcriptome reveals the lncRNA/circRNA regulation in tissue development. Rice 14:14. doi: 10.1186/s12284-021-00455-2

PubMed Abstract | CrossRef Full Text | Google Scholar