Next Article in Journal
Quantification of Phytochemicals in Cephalotaxus harringtonia: Insights into Plant Tissue-Specific Allocation
Previous Article in Journal
Microencapsulation of Bacillus megaterium in Humic Acid-Supplied Alginate Beads Enhances Tomato Growth and Suppresses the Root-Knot Nematode Meloidogyne javanica Under Greenhouse Conditions
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Horticulturae
Volume 10
Issue 12
10.3390/horticulturae10121285
horticulturae-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
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Transcriptome Profiling Implicates Non-Coding RNAs Involved in Flowering and Floral Organ Development in Water Lily

by
Hongbin Zhang
1,2,
Chengjun Zhao
1,2 and
Fei Chen
1,2,*
1
National Key Laboratory for Tropical Crop Breeding, School of Breeding and Multiplication, Sanya Institute of Breeding and Multiplication, Hainan University, Sanya 572025, China
2
College of Tropical Agriculture and Forestry, Hainan University, Sanya 572024, China
*
Author to whom correspondence should be addressed.
Horticulturae 2024, 10(12), 1285; https://doi.org/10.3390/horticulturae10121285
Submission received: 25 October 2024 / Revised: 27 November 2024 / Accepted: 29 November 2024 / Published: 2 December 2024
(This article belongs to the Section Floriculture, Nursery and Landscape, and Turf)
Browse Figures
Versions Notes

Abstract

:
In this study, we performed small RNA and whole-transcriptome sequencing of different tissues of Nymphaea minuta to systematically investigate the roles and regulatory mechanisms of miRNA, lncRNA, and circRNA in the regulation of flowering-related target genes. Fifteen samples were sequenced using the Illumina platform, with strict data quality control to ensure the reliability of the analysis. By applying multiple bioinformatics tools, miRNA, lncRNA, and circRNA were comprehensively identified, annotated, and functionally analyzed, with a focus on screening non-coding RNAs closely related to the flowering process. The results showed significant differential expression of these miRNAs and lncRNAs across different tissues, which influenced the expression of flowering-related genes through specific regulatory networks. The constructed gene co-expression network further revealed the central roles of these non-coding RNAs in flowering regulation. This study provides new insights into the flowering regulatory mechanisms of N. minuta, highlights the potential of this species for studying aquatic plant flowering mechanisms, and provides an important theoretical basis for gene function research in aquatic plants.

1. Introduction

The growth and development of plants, especially the flowering process, are regulated at multiple levels, with non-coding RNAs (such as miRNA, lncRNA, and circRNA) playing important roles as regulatory molecules widely involved in precise control of plant gene expression [1]. In flowering regulation, miRNAs and lncRNAs have significant effects on flowering time, floral organ formation, and flower maturation [2]. In recent years, with the rapid development of high-throughput sequencing technologies, there has been a deeper understanding of the biological functions of non-coding RNAs, particularly their roles in the regulatory networks of plant flowering [2,3]. Studies have shown that miRNAs, lncRNAs, and circRNAs not only play key roles in the overall development of plants but also significantly influence flowering behavior under different environmental conditions [3,4].
The water lily is an early flowering plant, considered a living fossil for studying floral development [5]. N. minuta is a small and fast-growing model water lily with important ecological and ornamental values, and it is characterized by its small flowers and strong reproductive ability, making it an ideal candidate for studying flowering mechanisms in aquatic plants [6]. However, the roles of small RNAs and long non-coding RNAs in N. minuta remain largely unexplored, representing a gap in current research on aquatic plant molecular biology [7]. The regulatory functions of miRNAs and lncRNAs, particularly in the context of flowering and floral organ development, have not yet been systematically studied in N. minuta, despite their potential to provide valuable insights into the regulatory processes of aquatic plant flowering.
To elucidate the regulatory roles of non-coding RNAs in different tissues of N. minuta during flowering, we conducted small RNA and whole-transcriptome sequencing of different tissues, including flowers, leaves, fruits, and seeds [8]. We systematically analyzed the expression characteristics of miRNA, lncRNA, and circRNA, and their regulatory mechanisms during flowering [9,10]. Through comprehensive bioinformatics analysis, we aimed to identify flowering-related target genes and characterize their potential regulatory networks involving non-coding RNAs. The core objective of this study was to determine the specific regulatory roles of miRNAs, lncRNAs, and circRNAs in the flowering process of N. minuta, thereby addressing the current research gap and providing a theoretical basis for the flowering gene regulatory network of this species. This study also offers a scientific reference for exploring flowering mechanisms in other aquatic plants [11]. The results of this study will contribute to a deeper understanding of the functions and roles of non-coding RNAs in plant flowering processes, particularly in the complexity of flowering regulation under specific environmental conditions, providing new directions and insights for basic research and applications in plant science [12].

2. Materials and Methods

2.1. Sample Selection and Sequencing

In this study, 15 samples, including three repeats of flower, fruit, seed, leaf, and root, of N. minuta were selected for small RNA and whole-transcriptome sequencing analysis. Small RNA sequencing was performed using the Illumina platform in Hiseq Single-End mode, and the data were stored in FASTQ format. Whole-transcriptome sequencing was conducted using the Illumina HiSeq platform with paired-end sequencing (2 × 150 bp).

2.2. Data Preprocessing

For both miRNA and whole-transcriptome sequencing data, quality control and filtering were necessary to ensure the accuracy and reliability of subsequent analyses. The Fastp tool was used to remove adapter-containing reads and filter data based on average base quality. A sliding window of five bases was used to trim reads with an average quality score below Q20, retaining high-quality sequences for subsequent analyses. The filtered sequences were then evaluated for length distribution, ranging from 18 to 36 nt. Duplicate sequences were removed, and their respective abundances and length distributions were recorded. The unique reads after deduplication were used for annotation and analysis [13].

2.3. Genomic Alignment

Both miRNA and whole-transcriptome data were aligned to the reference genome of the closely related species Nymphaea colorata (Genome GCF_008831285.2_ASM883128v2) [14]. Alignment was performed using software, such as Hisat2 (version 2.2.1) [15] and Bowtie (version 1.2.3) [16], with parameters set to allow a maximum of two mismatches to ensure the accuracy of alignment results and high retention of quality data.

2.4. miRNA Annotation and Expression Analysis

In miRNA annotation and analysis, the unique reads were classified and annotated using the miRBase database (v22) [17], the ncRNA database [18], the repeat database of N. minuta, intron sequences, and the mRNA sequences of N. colorata. Conserved miRNAs were screened, allowing for up to two mismatches for closely related species. For sequences that could not be annotated in the miRBase database, novel miRNAs were predicted using miRDeep-P2 software (v1.1.4) [19,20], which determined their locations and structural features, and RNAfold was used to draw secondary structure diagrams. The DESeq2 package (v1.46) [21] was used to calculate miRNA expression levels for each sample, with read counts normalized to TPM (Transcripts Per Million) for subsequent differential expression and clustering analyses. The psRobot software (v1.2) [22] was used for target gene prediction of differentially expressed miRNAs, and functional enrichment analysis was performed based on the Gene Ontology (GO) [23] and KEGG databases [24].

2.5. Whole-Transcriptome Analysis

In the whole-transcriptome analysis, gene expression levels were quantified using the RSEM software (v1.3.3) [25] based on clean reads, and expression was normalized using the FPKM (Fragments Per Kilobase Million) method. The StringTie software was used to assemble aligned reads and annotate transcripts. CNCI (v2) [26], PLEK(v2) [27], and CPC2 (v1) [28] were used to predict the non-coding potential of transcripts, and high-confidence lncRNAs were screened [29]. Unaligned reads were re-analyzed using find_circ [30] to identify and annotate circRNAs [31].

2.6. Differential Expression Analysis and Functional Enrichment Analysis

Differential expression analysis of miRNA, lncRNA, and circRNA across different tissues was conducted using DESeq2, with the screening criteria set as |log2 Fold Change| > 1 and p-value < 0.05. Cluster analysis was performed on differentially expressed genes using the Pheatmap package to evaluate expression patterns between samples. GO and KEGG enrichment analyses were performed on target genes of differentially expressed mRNAs and lncRNAs. GO enrichment analysis was divided into three categories: biological process (BP), cellular component (CC), and molecular function (MF). KEGG analysis was used to identify significantly enriched metabolic pathways. To further explore gene set-level functional enrichment, GSEA (Gene Set Enrichment Analysis) [32] was used for functional enrichment evaluation at the gene set level, revealing significant regulatory pathways in different tissues.

3. Results

3.1. Data Quality and miRNA Annotation Results

In this study, approximately 1.7 billion raw reads were obtained from 15 N. minuta samples (Figure 1a,b). After strict quality control, a total of 1.695 billion clean reads were retained, with an average Q30 value above 96%, indicating high data quality (Figure 1c,d). Annotation using the Rfam database identified four known classes of non-coding RNA (ncRNA), including rRNA, tRNA, snRNA, and snoRNA. After alignment and filtering, 3437 conserved miRNAs were screened from the 15 samples, and 696 unique conserved miRNAs were identified after deduplication (Figure 1e–g and Table 1). Some of these miRNAs showed significant differential expression across different tissues.

3.2. Novel miRNA Prediction and Functional Analysis

A total of 31 novel miRNAs were predicted using the miRDeep-P2 software, exhibiting significant differential expression among different tissues, particularly in flowers, leaves, and fruits. Target gene prediction of these miRNAs revealed their association with several important biological pathways, including plant hormone signaling, cell cycle regulation, and photosynthesis. Prior to the differential expression analysis, the correlation of microRNA expression levels between samples was examined, revealing a strong correlation between flowers, fruits, and seeds (Figure 2c). Differentially expressed conserved microRNAs were selected based on fold change and the statistical significance of expression differences (Figure 2d). Since genes of the same category may have similar functions, a bidirectional clustering analysis was performed on the genes and samples based on the differentially expressed genes across all groups, grouping genes with similar or close expression patterns together (Figure 2e,f). GO enrichment analysis showed that the target genes of these miRNAs were mainly involved in processes, such as photosynthesis, hormone signal transduction, and cell proliferation [33] (Figure 3a,b). KEGG analysis indicated that the target genes of these miRNAs were enriched in plant hormone signaling pathways and nitrogen metabolism pathways, highlighting their importance in plant growth and development [34] (Figure 3c,d). Due to the short length of microRNAs, only the total reads were used for inter-sample normalization. Overall gene expression patterns were examined using CPM density distribution plots, showing that moderately expressed miRNAs accounted for the majority, while low and high-expression miRNAs represented a small proportion [35] (Figure 2a,b).

3.3. Whole-Transcriptome Expression and Differential Analysis

For whole-transcriptome data, mRNA expression levels were quantified and standardized using FPKM values (Figure 4a). Principal component analysis (PCA) results showed distinct tissue-specific expression profiles for flowers, leaves, fruits, and seeds, indicating clear functional differences between these tissues (Figure 4b). A total of 1981 high-confidence lncRNAs were identified (Figure 4d,e; additional data available online). Prior to the differential expression analysis, the correlation of gene expression levels between samples was examined (Figure 4c). Differential expression analysis revealed that some lncRNAs were significantly upregulated in leaves, with their target genes mainly enriched in flowering-related biological processes, such as floral organ formation, petal differentiation, carbohydrate metabolism in floral structures, and circadian rhythm regulation associated with flowering [36]. GO and KEGG enrichment analyses indicated that the target genes of these lncRNAs were involved in metabolic pathways, starch and sucrose metabolism, glycolysis, and flowering-related circadian rhythm regulation [37] (Figure 4f,g).

3.4. circRNA Identification and Functional Prediction

A total of 1520 circRNAs were identified using the find_circ tool, most of which originated from intronic regions of the genome (Figure 5a–c; additional data available online). Functional enrichment analysis suggested that circRNAs might act as molecular sponges by competitively binding to miRNAs, thereby regulating downstream gene expression. This finding provides new insights into the role of circRNAs in plant gene expression regulation.

3.5. Functional Enrichment and Co-Expression Network Analysis of Differentially Expressed Genes

After performing functional enrichment analysis on differentially expressed miRNAs, lncRNAs, and circRNAs, GO analysis results indicated that these differentially expressed genes were mainly involved in biological processes related to flowering regulation, cell structure, and molecular transport [38]. KEGG analysis showed that these genes were significantly enriched in secondary metabolite synthesis pathways and flowering-related biological processes. Based on these data, a gene co-expression network, including mRNAs, lncRNAs, and circRNAs, was constructed. Analysis of the co-expression network suggested that lncRNAs and circRNAs may play important regulatory roles in floral organ development and flower maturation. From the enrichment analysis, we identified nine miRNAs and nine lncRNAs that regulate flowering-related target genes, which are involved in flowering regulation and floral organ development (Figure 6a,c). Among them, MSTRG.23579 was highly expressed in the regulation of flowering-related target genes (Figure 6b). These non-coding RNAs were closely associated with mRNAs in several key metabolic pathways, collectively participating in the regulation of floral organ formation and the flowering process.
Through comprehensive analysis of miRNA and whole-transcriptome sequencing data, this study revealed the gene expression patterns and underlying regulatory mechanisms in different tissues of N. minuta. The results indicated that miRNAs, lncRNAs, and circRNAs play significant regulatory roles in plant growth and development, particularly in leaf function and fruit maturation [39]. These findings provide a rich data foundation and theoretical basis for further research on plant developmental mechanisms, and future studies can build upon this work to explore the specific roles of non-coding RNAs in plant physiological functions and environmental adaptability.

4. Discussion

This study comprehensively identified and analyzed miRNAs, lncRNAs, and circRNAs in different tissues of N. minuta, revealing their important regulatory roles in flower development. First, the significant differential expression of miRNAs in different tissues suggests that they may have tissue-specific regulatory functions, particularly in the development of floral organs. Functional enrichment analysis of miRNA target genes showed that these genes were mainly involved in biological processes such as flowering regulation, hormone signal transduction, and cell cycle regulation, further supporting the critical role of miRNAs in regulating plant growth and stress responses.
lncRNAs and circRNAs were also shown to have important regulatory roles in this study. Differential expression analysis of lncRNAs revealed their significant tissue specificity during floral organ formation and maturation [40]. Specifically, some newly identified lncRNAs had target genes enriched in functional pathways related to floral organ development, petal formation, carbohydrate metabolism, and pyruvate metabolism, suggesting that these lncRNAs may play important regulatory roles in flower development. In addition, functional enrichment results for circRNAs indicated that they might act as “sponges” by competitively binding to miRNAs, thereby indirectly modulating mRNA expression. This competitive endogenous RNA (ceRNA) regulatory mechanism provides new insights into the complexity of gene expression regulation in plants, particularly in flower development [41].
Co-expression network analysis based on whole-transcriptome data further revealed the coordinated regulatory roles of miRNAs, lncRNAs, and circRNAs in the flower development process of plants. The constructed gene co-expression network showed that there are complex interactions between these non-coding RNAs and mRNAs, especially during flower formation and maturation, where multiple key metabolic and signaling pathways are jointly regulated by these non-coding RNAs. These results indicate that non-coding RNAs not only have regulatory roles at the individual gene level but also play important roles in higher-level gene networks, thereby achieving precise regulation of the complex biological processes involved in flower development.
Although this study revealed the expression patterns and potential regulatory mechanisms of miRNAs, lncRNAs, and circRNAs in N. minuta, several issues require further investigation. First, while a large number of non-coding RNAs and their target genes were predicted through bioinformatics analysis, experimental validation is lacking. Future research should incorporate experimental methods (such as qPCR and RNA interference) to verify the functions of these non-coding RNAs, further confirming their specific roles in plant growth and development. In addition, this study was limited to the analysis of specific tissues, and future studies could introduce more environmental factors and time points to explore the dynamic changes in non-coding RNAs under environmental stress and at different developmental stages.
In conclusion, this study used high-throughput sequencing and comprehensive bioinformatics analysis to reveal the expression patterns of miRNAs, lncRNAs, and circRNAs in N. minuta and their important regulatory roles in plant growth and development. The diversity and complexity of these non-coding RNAs in gene expression regulation provide new perspectives for understanding plant growth mechanisms and offer an important theoretical basis and data support for future gene function research in aquatic plants.

Author Contributions

F.C. designed and led this study. H.Z. and C.Z. performed the whole analyses and wrote the manuscript. All the authors approved the final manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (32172614), Hainan Province Science and Technology Special Fund (ZDYF2023XDNY050), Hainan Provincial Natural Science Foundation of China (324RC452), and the Project of National Key Laboratory for Tropical Crop Breeding (NO. NKLTCB202337).

Data Availability Statement

Raw reads have been submitted to the National Genomics Data Center, project ID:PRJCA031249 for miRNA and PRJCA031250 for whole-genome transcriptome. The miRNAs, lncRNAs, and circRNAs identified in this study are available in our database WLP, accessible at https://bioinformatics.hainanu.edu.cn/waterlily/ (accessed on 1 December 2024) under the Transcriptomics section.

Acknowledgments

We would like to express our sincere gratitude to the National Natural Science Foundation of China (32172614), Hainan Province Science and Technology Special Fund (ZDYF2023XDNY050), Hainan Provincial Natural Science Foundation of China (324RC452), and the Project of National Key Laboratory for Tropical Crop Breeding (NO. NKLTCB202337) for their generous financial support. Additionally, we would like to thank the editor and anonymous reviewers for their valuable feedback and insightful suggestions, which greatly improved the quality of this manuscript.

Conflicts of Interest

The authors declare that there are no conflicts of interest.

References

  1. Zhao, Z.; Zang, S.; Zou, W.; Pan, Y.B.; Yao, W.; You, C.; Que, Y. Long Non-Coding RNAs: New Players in Plants. Int. J. Mol. Sci. 2022, 23, 9301. [Google Scholar] [CrossRef] [PubMed]
  2. Jeena, G.S.; Singh, N.; Shikha; Shukla, R.K. An insight into microRNA biogenesis and its regulatory role in plant secondary metabolism. Plant Cell Rep. 2022, 41, 1651–1671. [Google Scholar] [CrossRef] [PubMed]
  3. Zhang, Q.; Zhao, Y.Q.; Gao, X.; Jia, G.X. Analysis of miRNA-mediated regulation of flowering induction in Lilium x formolongi. BMC Plant Biol. 2021, 21, 190. [Google Scholar] [CrossRef] [PubMed]
  4. Xu, J.; Wang, Q.; Tang, X.; Feng, X.; Zhang, X.; Liu, T.; Wu, F.; Wang, Q.; Feng, X.; Tang, Q.; et al. Drought-induced circular RNAs in maize roots: Separating signal from noise. Plant Physiol. 2024, 196, 352–367. [Google Scholar] [CrossRef] [PubMed]
  5. Xiong, X.; Zhang, J.; Yang, Y.; Chen, Y.; Su, Q.; Zhao, Y.; Wang, J.; Xia, Z.; Wang, L.; Zhang, L.; et al. Water lily research: Past, present, and future. Trop. Plants 2023, 2, 1. [Google Scholar] [CrossRef]
  6. Chen, F.; Liu, X.; Yu, C.; Chen, Y.; Tang, H.; Zhang, L. Water lilies as emerging models for Darwin’s abominable mystery. Hortic. Res. 2017, 4, 17051. [Google Scholar] [CrossRef]
  7. Zhang, Y.; Rahmani, R.S.; Yang, X.; Chen, J.; Shi, T. Integrative expression network analysis of microRNA and gene isoforms in sacred lotus. BMC Genom. 2020, 21, 429. [Google Scholar] [CrossRef]
  8. Yan, J.; Gu, Y.; Jia, X.; Kang, W.; Pan, S.; Tang, X.; Chen, X.; Tang, G. Effective small RNA destruction by the expression of a short tandem target mimic in Arabidopsis. Plant Cell 2012, 24, 415–427. [Google Scholar] [CrossRef]
  9. Palos, K.; Nelson Dittrich, A.C.; Yu, L.; Brock, J.R.; Railey, C.E.; Wu, H.L.; Sokolowska, E.; Skirycz, A.; Hsu, P.Y.; Gregory, B.D.; et al. Identification and functional annotation of long intergenic non-coding RNAs in Brassicaceae. Plant Cell 2022, 34, 3233–3260. [Google Scholar] [CrossRef]
  10. Zhou, Y.; Myat, A.A.; Liang, C.; Meng, Z.; Guo, S.; Wei, Y.; Sun, G.; Wang, Y.; Zhang, R. Insights Into MicroRNA-Mediated Regulation of Flowering Time in Cotton Through Small RNA Sequencing. Front. Plant Sci. 2022, 13, 761244. [Google Scholar] [CrossRef]
  11. Bernardi, Y.; Ponso, M.A.; Belen, F.; Vegetti, A.C.; Dotto, M.C. MicroRNA miR394 regulates flowering time in Arabidopsis thaliana. Plant Cell Rep. 2022, 41, 1375–1388. [Google Scholar] [CrossRef] [PubMed]
  12. Zhang, L.; Lin, T.; Zhu, G.; Wu, B.; Zhang, C.; Zhu, H. LncRNAs exert indispensable roles in orchestrating the interaction among diverse noncoding RNAs and enrich the regulatory network of plant growth and its adaptive environmental stress response. Hortic. Res. 2023, 10, uhad234. [Google Scholar] [CrossRef] [PubMed]
  13. Zhan, J.; Meyers, B.C. Plant Small RNAs: Their Biogenesis, Regulatory Roles, and Functions. Annu. Rev. Plant Biol. 2023, 74, 21–51. [Google Scholar] [CrossRef] [PubMed]
  14. Zhang, L.; Chen, F.; Zhang, X.; Li, Z.; Zhao, Y.; Lohaus, R.; Chang, X.; Dong, W.; Ho, S.Y.W.; Liu, X.; et al. The water lily genome and the early evolution of flowering plants. Nature 2020, 577, 79–84. [Google Scholar] [CrossRef]
  15. Kim, D.; Langmead, B.; Salzberg, S.L. HISAT: A fast spliced aligner with low memory requirements. Nat. Methods 2015, 12, 357–360. [Google Scholar] [CrossRef]
  16. Langmead, B.; Salzberg, S.L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 2012, 9, 357–359. [Google Scholar] [CrossRef]
  17. Kozomara, A.; Birgaoanu, M.; Griffiths-Jones, S. miRBase: From microRNA sequences to function. Nucleic Acids Res. 2019, 47, D155–D162. [Google Scholar] [CrossRef]
  18. Zheng, Y.; Luo, H.; Teng, X.; Hao, X.; Yan, X.; Tang, Y.; Zhang, W.; Wang, Y.; Zhang, P.; Li, Y.; et al. NPInter v5.0: ncRNA interaction database in a new era. Nucleic Acids Res. 2023, 51, D232–D239. [Google Scholar] [CrossRef]
  19. Wang, Y.; Kuang, Z.; Li, L.; Yang, X. A Bioinformatics Pipeline to Accurately and Efficiently Analyze the MicroRNA Transcriptomes in Plants. J. Vis. Exp. 2020, 155, e59864. [Google Scholar] [CrossRef]
  20. Kuang, Z.; Wang, Y.; Li, L.; Yang, X. miRDeep-P2: Accurate and fast analysis of the microRNA transcriptome in plants. Bioinformatics 2019, 35, 2521–2522. [Google Scholar] [CrossRef]
  21. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550. [Google Scholar] [CrossRef] [PubMed]
  22. Wu, H.J.; Ma, Y.K.; Chen, T.; Wang, M.; Wang, X.J. PsRobot: A web-based plant small RNA meta-analysis toolbox. Nucleic Acids Res. 2012, 40, W22–W28. [Google Scholar] [CrossRef] [PubMed]
  23. Gene Ontology Consortium: Going forward. Nucleic Acids Res. 2015, 43, D1049–D1056. [CrossRef]
  24. Kanehisa, M.; Goto, S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000, 28, 27–30. [Google Scholar] [CrossRef]
  25. Li, B.; Dewey, C.N. RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 2011, 12, 323. [Google Scholar] [CrossRef]
  26. Sun, L.; Luo, H.; Bu, D.; Zhao, G.; Yu, K.; Zhang, C.; Liu, Y.; Chen, R.; Zhao, Y. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013, 41, e166. [Google Scholar] [CrossRef]
  27. Li, A.; Zhang, J.; Zhou, Z. PLEK: A tool for predicting long non-coding RNAs and messenger RNAs based on an improved k-mer scheme. BMC Bioinform. 2014, 15, 311. [Google Scholar] [CrossRef]
  28. Kang, Y.J.; Yang, D.C.; Kong, L.; Hou, M.; Meng, Y.Q.; Wei, L.; Gao, G. CPC2: A fast and accurate coding potential calculator based on sequence intrinsic features. Nucleic Acids Res. 2017, 45, W12–W16. [Google Scholar] [CrossRef]
  29. Tian, X.C.; Chen, Z.Y.; Nie, S.; Shi, T.L.; Yan, X.M.; Bao, Y.T.; Li, Z.C.; Ma, H.Y.; Jia, K.H.; Zhao, W.; et al. Plant-LncPipe: A computational pipeline providing significant improvement in plant lncRNA identification. Hortic. Res. 2024, 11, uhae041. [Google Scholar] [CrossRef]
  30. Memczak, S.; Jens, M.; Elefsinioti, A.; Torti, F.; Krueger, J.; Rybak, A.; Maier, L.; Mackowiak, S.D.; Gregersen, L.H.; Munschauer, M.; et al. Circular RNAs are a large class of animal RNAs with regulatory potency. Nature 2013, 495, 333–338. [Google Scholar] [CrossRef]
  31. Chen, L.; Wang, C.; Sun, H.; Wang, J.; Liang, Y.; Wang, Y.; Wong, G. The bioinformatics toolbox for circRNA discovery and analysis. Brief. Bioinform. 2021, 22, 1706–1728. [Google Scholar] [CrossRef] [PubMed]
  32. Powers, R.K.; Goodspeed, A.; Pielke-Lombardo, H.; Tan, A.C.; Costello, J.C. GSEA-InContext: Identifying novel and common patterns in expression experiments. Bioinformatics 2018, 34, i555–i564. [Google Scholar] [CrossRef] [PubMed]
  33. Zhang, S.; Wu, Y.; Huang, X.; Wu, W.; Lyu, L.; Li, W. Research progress about microRNAs involved in plant secondary metabolism. Int. J. Biol. Macromol. 2022, 216, 820–829. [Google Scholar] [CrossRef] [PubMed]
  34. Xing, L.; Zhang, D.; Zhao, C.; Li, Y.; Ma, J.; An, N.; Han, M. Shoot bending promotes flower bud formation by miRNA-mediated regulation in apple (Malus domestica Borkh.). Plant Biotechnol. J. 2016, 14, 749–770. [Google Scholar] [CrossRef]
  35. Franco-Zorrilla, J.M.; Valli, A.; Todesco, M.; Mateos, I.; Puga, M.I.; Rubio-Somoza, I.; Leyva, A.; Weigel, D.; Garcia, J.A.; Paz-Ares, J. Target mimicry provides a new mechanism for regulation of microRNA activity. Nat. Genet. 2007, 39, 1033–1037. [Google Scholar] [CrossRef]
  36. Xu, S.; Dong, Q.; Deng, M.; Lin, D.; Xiao, J.; Cheng, P.; Xing, L.; Niu, Y.; Gao, C.; Zhang, W.; et al. The vernalization-induced long non-coding RNA VAS functions with the transcription factor TaRF2b to promote TaVRN1 expression for flowering in hexaploid wheat. Mol. Plant 2021, 14, 1525–1538. [Google Scholar] [CrossRef]
  37. Yadav, A.; Mathan, J.; Dubey, A.K.; Singh, A. The Emerging Role of Non-Coding RNAs (ncRNAs) in Plant Growth, Development, and Stress Response Signaling. Noncoding RNA 2024, 10, 13. [Google Scholar] [CrossRef]
  38. Diener, C.; Keller, A.; Meese, E. The miRNA-target interactions: An underestimated intricacy. Nucleic Acids Res. 2024, 52, 1544–1557. [Google Scholar] [CrossRef]
  39. Jin, Y.; Ivanov, M.; Dittrich, A.N.; Nelson, A.D.; Marquardt, S. LncRNA FLAIL affects alternative splicing and represses flowering in Arabidopsis. EMBO J. 2023, 42, e110921. [Google Scholar] [CrossRef]
  40. Xue, S.T.; Zheng, B.; Cao, S.Q.; Ding, J.C.; Hu, G.S.; Liu, W.; Chen, C. Long non-coding RNA LINC00680 functions as a ceRNA to promote esophageal squamous cell carcinoma progression through the miR-423-5p/PAK6 axis. Mol. Cancer 2022, 21, 69. [Google Scholar] [CrossRef]
  41. Teotia, S.; Tang, G. To bloom or not to bloom: Role of microRNAs in plant flowering. Mol. Plant 2015, 8, 359–377. [Google Scholar] [CrossRef]
Horticulturae 10 01285 g001
Figure 1. Statistics of microRNA from N. minuta. (a) Length distribution of each sequence across 15 samples after removing duplicates. (b) Number of shared and unique sequences among different groups after removing duplicates. (c) First-base preference of a single sample (flower1) after removing duplicates. The x-axis represents different miRNA lengths, and the y-axis indicates the percentage frequency of different bases at the first position of miRNAs with various lengths. (d) Base preference at each position of a single sample (flower1) after removing duplicates. The x-axis represents the positions of miRNAs, and the y-axis indicates the percentage frequency of different bases at each position of the miRNAs. (e) Nucleotide preference at the first position of 696 conserved miRNAs. (f) Nucleotide preference at each position of 696 conserved miRNAs. (g) Length distribution of 696 conserved-miRNA sequences.
Figure 1. Statistics of microRNA from N. minuta. (a) Length distribution of each sequence across 15 samples after removing duplicates. (b) Number of shared and unique sequences among different groups after removing duplicates. (c) First-base preference of a single sample (flower1) after removing duplicates. The x-axis represents different miRNA lengths, and the y-axis indicates the percentage frequency of different bases at the first position of miRNAs with various lengths. (d) Base preference at each position of a single sample (flower1) after removing duplicates. The x-axis represents the positions of miRNAs, and the y-axis indicates the percentage frequency of different bases at each position of the miRNAs. (e) Nucleotide preference at the first position of 696 conserved miRNAs. (f) Nucleotide preference at each position of 696 conserved miRNAs. (g) Length distribution of 696 conserved-miRNA sequences.
Horticulturae 10 01285 g001
Horticulturae 10 01285 g002
Figure 2. microRNA Expression Analysis. (a) CPM density plot. CPM is calculated as CPM = C/N × 1,000,000, where C is the number of reads mapped to a gene, and N is the total number of reads mapped to all genes. This plot examines the gene expression patterns of the samples. (b) CPM violin plot. The horizontal line within the box represents the median, the box edges represent the 75th percentile, and the upper/lower whiskers indicate the 90th percentile. The outer shape reflects kernel density estimation. (c) Correlation analysis of 15 sample groups. The left and top sections show sample clustering, while the right and bottom sections display sample names. The color intensity of the squares indicates the degree of correlation between samples. Generally, a correlation coefficient between 0.8 and 1 is considered to indicate a very strong correlation. (d) Volcano and MA plots of differentially expressed genes between leaf and flower. There are 44 low-expression genes (blue), 44 high-expression genes (red), and 395 non-significant genes (gray). In the volcano plot, the x-axis represents the log2 fold-change, and the y-axis shows the −log10 p-value. In the MA plot, the x-axis is the sum of gene expression (log2(A) + log2(B)), where A and B are the expression levels in two samples, and the y-axis is the expression difference (log2(A) − log2(B)). (e) Heatmap of differentially expressed miRNAs across tissues. Genes with similar or identical expression patterns are clustered to identify unknown functions of known genes or functions of novel genes, as co-expressed genes are likely functionally related. (f) Clustering analysis based on hierarchical clustering results. Differentially expressed genes are grouped into clusters according to expression patterns. Each gray line represents the expression profile of a gene within a cluster, while the blue line shows the average expression level of all genes in the cluster across samples.
Figure 2. microRNA Expression Analysis. (a) CPM density plot. CPM is calculated as CPM = C/N × 1,000,000, where C is the number of reads mapped to a gene, and N is the total number of reads mapped to all genes. This plot examines the gene expression patterns of the samples. (b) CPM violin plot. The horizontal line within the box represents the median, the box edges represent the 75th percentile, and the upper/lower whiskers indicate the 90th percentile. The outer shape reflects kernel density estimation. (c) Correlation analysis of 15 sample groups. The left and top sections show sample clustering, while the right and bottom sections display sample names. The color intensity of the squares indicates the degree of correlation between samples. Generally, a correlation coefficient between 0.8 and 1 is considered to indicate a very strong correlation. (d) Volcano and MA plots of differentially expressed genes between leaf and flower. There are 44 low-expression genes (blue), 44 high-expression genes (red), and 395 non-significant genes (gray). In the volcano plot, the x-axis represents the log2 fold-change, and the y-axis shows the −log10 p-value. In the MA plot, the x-axis is the sum of gene expression (log2(A) + log2(B)), where A and B are the expression levels in two samples, and the y-axis is the expression difference (log2(A) − log2(B)). (e) Heatmap of differentially expressed miRNAs across tissues. Genes with similar or identical expression patterns are clustered to identify unknown functions of known genes or functions of novel genes, as co-expressed genes are likely functionally related. (f) Clustering analysis based on hierarchical clustering results. Differentially expressed genes are grouped into clusters according to expression patterns. Each gray line represents the expression profile of a gene within a cluster, while the blue line shows the average expression level of all genes in the cluster across samples.
Horticulturae 10 01285 g002
Horticulturae 10 01285 g003
Figure 3. GO and KEGG Enrichment Analysis of Differential microRNA Target Genes in Different Parts of Water Lily (Leaf and Flower Displayed). (a) Bubble plot of GO enrichment analysis for target genes of differential microRNAs, showing the top 20 GO terms with the smallest FDR values (most significant enrichment). (b) Bar plot of GO enrichment analysis for target genes of differential microRNAs, categorized by MF, BP, and CC, displaying the top 10 GO terms with the smallest p-values (most significant enrichment) in each category. (c) Bubble plot of KEGG enrichment analysis for target genes of differential microRNAs, showing the top 20 KEGG pathways with the smallest FDR values (most significant enrichment). Two pathways are enriched (red) in the figure. Enrichment of both pathways typically indicates that the cell is in a state requiring heightened protein synthesis and gene expression regulation. (d) Bar plot of KEGG enrichment analysis for target genes of differential microRNAs, displaying the top 30 pathways with the smallest p-values (most significant enrichment).
Figure 3. GO and KEGG Enrichment Analysis of Differential microRNA Target Genes in Different Parts of Water Lily (Leaf and Flower Displayed). (a) Bubble plot of GO enrichment analysis for target genes of differential microRNAs, showing the top 20 GO terms with the smallest FDR values (most significant enrichment). (b) Bar plot of GO enrichment analysis for target genes of differential microRNAs, categorized by MF, BP, and CC, displaying the top 10 GO terms with the smallest p-values (most significant enrichment) in each category. (c) Bubble plot of KEGG enrichment analysis for target genes of differential microRNAs, showing the top 20 KEGG pathways with the smallest FDR values (most significant enrichment). Two pathways are enriched (red) in the figure. Enrichment of both pathways typically indicates that the cell is in a state requiring heightened protein synthesis and gene expression regulation. (d) Bar plot of KEGG enrichment analysis for target genes of differential microRNAs, displaying the top 30 pathways with the smallest p-values (most significant enrichment).
Horticulturae 10 01285 g003
Horticulturae 10 01285 g004
Figure 4. lncRNA Expression Levels and Differential Expression Analysis. (a) FPKM density distribution plot to examine the expression patterns of lncRNAs across 15 samples. The majority of lncRNAs exhibit moderate expression levels. The x-axis represents the log10FPKM values of genes, while the y-axis shows the density distribution of lncRNAs at corresponding expression levels. (b) Principal Component Analysis (PCA) of lncRNA expression levels across samples using DESeq software. PCA clusters similar samples close together, with shorter distances indicating higher similarity between samples. Different shapes represent different samples, and different colors correspond to different groups. (c) Pearson correlation coefficients between samples indicate the correlation of lncRNA expression levels. The closer the correlation coefficient is to 1, the higher the similarity between the expression patterns of the samples. (d) Comparison of the distribution of exon numbers between lncRNA and mRNA. The x-axis represents the number of exons in the sequences, and the y-axis represents the number of sequences that contain that number of exons. (e) Venn analysis using the VennDiagram package to compare non-coding potential transcripts predicted by PLEK, CNCI, CPC2, and Blast. Transcripts identified by all four tools as lacking coding potential are considered high-confidence lncRNAs and are used for further analysis. (f) GO enrichment analysis of lncRNA target genes, showing the top 20 GO terms with the smallest FDR values, indicating the most significantly enriched terms. The target genes are mainly associated with regions such as thylakoids and photosynthetic membranes. (g) KEGG enrichment analysis of lncRNA target genes, displaying the top 20 pathways with the smallest FDR values, representing the most significantly enriched pathways. These pathways are primarily involved in photosynthesis, pyruvate metabolism, and starch and sucrose metabolism.
Figure 4. lncRNA Expression Levels and Differential Expression Analysis. (a) FPKM density distribution plot to examine the expression patterns of lncRNAs across 15 samples. The majority of lncRNAs exhibit moderate expression levels. The x-axis represents the log10FPKM values of genes, while the y-axis shows the density distribution of lncRNAs at corresponding expression levels. (b) Principal Component Analysis (PCA) of lncRNA expression levels across samples using DESeq software. PCA clusters similar samples close together, with shorter distances indicating higher similarity between samples. Different shapes represent different samples, and different colors correspond to different groups. (c) Pearson correlation coefficients between samples indicate the correlation of lncRNA expression levels. The closer the correlation coefficient is to 1, the higher the similarity between the expression patterns of the samples. (d) Comparison of the distribution of exon numbers between lncRNA and mRNA. The x-axis represents the number of exons in the sequences, and the y-axis represents the number of sequences that contain that number of exons. (e) Venn analysis using the VennDiagram package to compare non-coding potential transcripts predicted by PLEK, CNCI, CPC2, and Blast. Transcripts identified by all four tools as lacking coding potential are considered high-confidence lncRNAs and are used for further analysis. (f) GO enrichment analysis of lncRNA target genes, showing the top 20 GO terms with the smallest FDR values, indicating the most significantly enriched terms. The target genes are mainly associated with regions such as thylakoids and photosynthetic membranes. (g) KEGG enrichment analysis of lncRNA target genes, displaying the top 20 pathways with the smallest FDR values, representing the most significantly enriched pathways. These pathways are primarily involved in photosynthesis, pyruvate metabolism, and starch and sucrose metabolism.
Horticulturae 10 01285 g004
Horticulturae 10 01285 g005
Figure 5. High-Confidence CircRNA Identification and Genomic Distribution. (a) Length distribution characteristics of CircRNAs across samples. CircRNAs were identified using the find_circ tool, with the following filtering criteria: anchor_overlap ≤ 2, allowing a maximum of 2 bp mismatches, and CircRNA length < 100 kb. For each read, at least one of the two anchor sequences must have a genome alignment score higher than 35. (b) Distribution of CircRNA types across samples. Annot_exon or one_exon: CircRNAs composed of a single exon (one_exon) or multiple exons (annot_exon). Exon_intron: CircRNAs containing one or more exons, but the breakpoint does not fall within an exon region. Intronic: CircRNAs located within introns without including any known exons. Antisense: CircRNAs found on the antisense strand of known genes. Intergenic: CircRNAs located between known genes. (c) Chromosomal distribution of CircRNA counts across samples, showing how many CircRNAs are mapped to each chromosome.
Figure 5. High-Confidence CircRNA Identification and Genomic Distribution. (a) Length distribution characteristics of CircRNAs across samples. CircRNAs were identified using the find_circ tool, with the following filtering criteria: anchor_overlap ≤ 2, allowing a maximum of 2 bp mismatches, and CircRNA length < 100 kb. For each read, at least one of the two anchor sequences must have a genome alignment score higher than 35. (b) Distribution of CircRNA types across samples. Annot_exon or one_exon: CircRNAs composed of a single exon (one_exon) or multiple exons (annot_exon). Exon_intron: CircRNAs containing one or more exons, but the breakpoint does not fall within an exon region. Intronic: CircRNAs located within introns without including any known exons. Antisense: CircRNAs found on the antisense strand of known genes. Intergenic: CircRNAs located between known genes. (c) Chromosomal distribution of CircRNA counts across samples, showing how many CircRNAs are mapped to each chromosome.
Horticulturae 10 01285 g005
Horticulturae 10 01285 g006
Figure 6. Functional enrichment analysis of target genes related to flowering regulation and floral organ development in N. minuta. (a) Functional enrichment of target genes involved in flowering regulation and floral organ development among miRNAs and lncRNAs. No flower-related target genes were identified for circRNAs. (b) Expression patterns of nine flowering-related miRNAs and nine lncRNAs across 15 sample groups. (c) Comparison between miRNAs and lncRNAs and their corresponding regulated target genes. The parts marked in red in the figure represent the highly expressed genes.
Figure 6. Functional enrichment analysis of target genes related to flowering regulation and floral organ development in N. minuta. (a) Functional enrichment of target genes involved in flowering regulation and floral organ development among miRNAs and lncRNAs. No flower-related target genes were identified for circRNAs. (b) Expression patterns of nine flowering-related miRNAs and nine lncRNAs across 15 sample groups. (c) Comparison between miRNAs and lncRNAs and their corresponding regulated target genes. The parts marked in red in the figure represent the highly expressed genes.
Horticulturae 10 01285 g006
Table 1. Number of microRNAs identified in 15 transcriptome samples.
Table 1. Number of microRNAs identified in 15 transcriptome samples.
SamplemiRNAUniqueTotal
flower129919041,071,228
flower22921873976,862
flower328218531,063,654
fruit12021130419,238
fruit2180939286,935
fruit3187891236,114
leaf127619952,590,514
leaf226619502,144,999
leaf326118402,104,628
root12091188434,006
root21981185442,464
root328517241,195,523
seed116952534,734
seed217456739,231
seed315751733,604
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

Zhang, H.; Zhao, C.; Chen, F. Transcriptome Profiling Implicates Non-Coding RNAs Involved in Flowering and Floral Organ Development in Water Lily. Horticulturae 2024, 10, 1285. https://doi.org/10.3390/horticulturae10121285

AMA Style

Zhang H, Zhao C, Chen F. Transcriptome Profiling Implicates Non-Coding RNAs Involved in Flowering and Floral Organ Development in Water Lily. Horticulturae. 2024; 10(12):1285. https://doi.org/10.3390/horticulturae10121285

Chicago/Turabian Style

Zhang, Hongbin, Chengjun Zhao, and Fei Chen. 2024. "Transcriptome Profiling Implicates Non-Coding RNAs Involved in Flowering and Floral Organ Development in Water Lily" Horticulturae 10, no. 12: 1285. https://doi.org/10.3390/horticulturae10121285

APA Style

Zhang, H., Zhao, C., & Chen, F. (2024). Transcriptome Profiling Implicates Non-Coding RNAs Involved in Flowering and Floral Organ Development in Water Lily. Horticulturae, 10(12), 1285. https://doi.org/10.3390/horticulturae10121285

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