Transcriptomic Insight into Viviparous Growth in Water Lily

Water lily is an important ornamental flower plant which is capable of viviparous plantlet development. But no study has been reported on the molecular basis of viviparity in water lily. Hence, we performed a comparative transcriptome study between viviparous water lily Nymphaea micrantha and a nonviviparous species Nymphaea colorata at four developmental stages. The higher expression of highly conserved AUX/IAA, ARF, GH3, and SAUR gene families in N. micrantha compared to N. colorata is predicted to have a major impact on the development and evolution of viviparity in water lily. Likewise, differential regulation of hormone signaling, brassinosteroid, photosynthesis, and energy-related pathways in the two species provide clues of their involvement in viviparity phenomenon. This study revealed the complex mechanism of viviparity trait in water lily. The transcriptomic signatures identified are important basis for future breeding and research of viviparity in water lily and other plant species.


Introduction
Water lilies (Nymphaea spp.) are ornamental flower plants distributed worldwide from the tropics to temperate regions. They are well known and widely cultivated for environment beautification as well as water purification [1]. They have rich flower colors, long flowering periods, strong adaptability, and stress resistance and are easy to grow. As a precious flower among aquatic flowers, water lily is more popular with enthusiast breeders and botanists because of its charming appearance and unique taxonomic status. Water lilies can be classified as hardy water lilies and tropical water lilies on the basis of the typical characteristics in different ecosystems [2]. The demand for water lilies with specific traits is rapidly increasing; hence, it is required to breed new water lily varieties and hybrids with excellent ornamental characteristics. The breeders are making efforts to develop new cultivars with artificial interspecific hybridization, mutation techniques, and polyploidy approaches. Among these tech-niques, the hybridization is the most widely adopted method [3]. However, the existence of incongruity barriers and the development of nonviable seeds reduce the breeding efficiency and speed of multiplication [2,3].
The breeding through contrasting parents is a conventional tool for cultivar development in plants. In some water lily species, their leaves along with normal photosynthesis also have a reproductive function known as viviparity ( Figure 1) [4]. Viviparous water lilies have the capacity to produce new plants that emerge while still attached to the parent plant known as plantlets. The slightly concave or smooth nose-like structure at the junction of the stem and leaf grows new plantlet, which can reproduce asexually in a large number of individuals in a short period of time [5].
There are a few tropical day blooming water lilies that produce plantlets from their pads. A few of the tropical night bloomers produce these plantlets from flowers [5]. The tropical day blooming variety Nymphaea micrantha has the ability to produce plantlets from its leaves, and it is often included in the breeding of all existing viviparous cultivars. The propagation of viviparous tropical water lilies can exponentially hasten their growth. The viviparous plantlets become mature in 20-30 days if cultivated in controlled conditions [6]. Recently, some studies on whole genome sequencing of water lilies have been reported [7,8]. Over the past few years, application of plant molecular biotechnological tools such as pollen-tube transgenics led to a lot of achievements in improving cold stress resistance and other traits in water lilies [9]. Nonetheless, the genetic mechanism and internal metabolism of water lily leaf viviparous phenomenon are unclear.
Currently, researchers mainly focus on evolution and taxonomic status, regulation of flower color, floral fragrance, flowering period, and reproduction of water lilies [10]. With the progress in genome sequencing tools, the use of transcriptome technology may provide a new way for mining the genes related to viviparous growth in water lilies [10]. In particular, transcriptome sequencing is a useful method for identifying novel transcripts and analyzing gene expression [11]. Through RNA-seq technology, a large amount of transcript data can be generated and manipulated to evaluate the gene expression, function, and related metabolic pathways. RNA-seq or transcriptome studies have extensively been used to evaluate many plant species [12,13] for various traits including early maturity [14] and seed dormancy [15]. The viviparity trait was investigated based on transcriptome analysis in mangrove plants [15,16] and Ophiorrhiza mungos L. [4], but limited transcriptome data exists in water lilies [10].
To understand the molecular mechanism of viviparous development in water lilies, Illumina paired-end sequencing analysis of the leaf tissues from N. micrantha with viviparous leaves and N. colorata with no viviparous leaves was conducted. It is not possible to make intraspecies comparison for this trait. We reported the differentially expressed genes and their functional enrichment between both species. This comprehensive analysis of the transcriptome substantially improved the understanding of the molecular mechanisms underlying viviparous growth in plant.

Plant Materials, Growth Conditions, Morphological
Identification, and RNA Extraction. The plants of two Nymphaea species N. micrantha (denoted with "X" in this study) and N. colorata (denoted with "L" in this study) (Figures 2(a) and 2(b)) were grown in a planting container of 100 cm diameter and 60 cm depth, with 30 cm soil thickness with standard natural growth conditions at Water Lily Resource Garden of the Floriculture Institute of Guangxi Academy of Agricultural Sciences, Guangxi, China, during the growth season of 2019. Six disease-free plants with robust growth for each genotype were selected for the transcriptome study. The morphological identification was performed when the leaves were fully unfolded and a clear viviparous mark in center could be observed. The leaf sampling was performed at four growth stages as (1) life size 1-2 cm long, (2) leaf size 4-6 cm, (3) when leaf tip connected to water surface, and (4) fully developed leaves. The fresh leaf samples with the leaf nose parts were collected at four developmental stages (Figures 2(c) and 2(d)). The leaves from six plants of the same genotype as technical replicates were mixed for each sample to make a homogenized sample, and a total of three biological repeats were obtained. The harvested leaves were immediately frozen in liquid nitrogen and stored at -80°C.
2.2. RNA Extraction and Preparation of Library. The total RNA was extracted by CTAB method [17] from leaf tissues of each sample using the three biological repeats. The extracted RNA was further evaluated for concentration (by NanoDrop 2000 microspectrophotometer), purity (by Agilent 2100 Bioanalyzer), and integrity in agarose gel. The step-by-step process as total RNA sample detection, mRNA enrichment with Oligo (dT) beads, fragmentation with   BioMed Research International fragmentation buffer, synthesis of double-stranded cDNA, purification and end repair, splice selection and PCR amplification library quality detection, and the computer-based sequencing was performed at Entrusted Jiugui Biotechnology Company Limited, Shanghai, China, and finally, the paired-end reads were generated.
2.3. Transcriptome Sequencing, Cleaning, and Assembly. The original image data files were obtained by high-throughput sequencing (Illumina HiSeq 4000) and were transformed into raw reads by base calling analysis. As per machine's sequencing strategy, 150 bp average read length was maintained. Raw data were processed using NGS QC Toolkit [18]. The raw reads with joint sequences, and/or less than 5 mass value, more than or equal to 50% proportion rate, more than or equal to 5% N-base (the base with undetermined information), containing Poly-A were filtered out to get the cleaned reads. Because the reference genome of N. micrantha is not yet available, we decided to a de novo assembly of the transcriptomes of the two species in this study. The Trinity v 2.6.6 program [19] was used for transcriptome assembly and to get the unigenes. The accuracy and effectiveness of the assembly results were ensured by estimating the N50 and exN50.

Expression Evaluation and Identification of Differentially
Expressed Genes. The number of reads count on each gene was obtained from each sample, and the gene expression level was estimated by the fragments per kilobase of transcript per million mapped reads (FPKM) method. FPKM value of each gene was calculated using cufflinks [20], and the read counts of each gene were obtained by htseq-count [21]. Differentially expressed genes (DEGs) were identified using the DESeq (with replicates) [22]. P value < 0.05 and log2 fold change > 1 for upregulated and fold change < −1 for downregulated DEGs were set as the threshold for significant differential expression. Principal component analysis of DEGs was performed to explore the gene expression pattern.

Functional
Annotation and Enrichment Analysis. The extracted unigenes were manipulated by Transcoder software v 4.1.0 to predict and translate the reading frames. Gene Ontology (GO) [23] enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway [24] enrichment analysis of DEGs were, respectively, performed using R based on the hypergeometric distribution. Blast2go [25] and Kaas software [26] (https://www.genome.jp/tools/kaas/) tools were used for GO and KEGG annotation, respectively, and Phyper function in R software was used for enrichment analysis. The gene expression was determined by comparing the sequenced reads with the unigene library in Bowtie [27].

Expression Validation by Quantitative Real-Time PCR.
Ten DEGs were further evaluated by qRT-PCR. The quantitative real-time PCR was performed using SYBR green mixture on an ABI 7500 real-time PCR detection system following the descriptions of Komivi et al. [28]. The Ubiquitin gene was used as internal control for normalization. The statistical approach developed by Livak and Schmittgen [29] was employed for statistical analysis. The expression data was further evaluated by Student's t-test for significance estimation. All primer sequences are listed in Additional  Table S4.

Results
3.1. Morphological Indication of Viviparity. The leaf tissues were evaluated at four developmental stages of the viviparous species N. micrantha denoted as X and nonviviparous species N. colorata denoted as L ( Figure 2). The four developmental stages can be defined as follows: stage (1) when the leaves are submerged and completely rolled with 1-2 cm length, stage (2) when the leaves are submerged and completely rolled with 4-6 cm length, stage (3) when the leaves are half rolled with the tip in contact with water     biological repeats. The high-throughput sequencing (Illumina HiSeq 4000 platform) data was generated and then transformed into the raw data by base calling analysis. A maximum 58.69 and 61.75 million raw reads were extracted for X and L, respectively (   Flavonoid biosynthesis 0.00000109 10 ko00196 Photosynthesis-antenna proteins 0.0000184 4 ko00565 Ether lipid metabolism 0.0000184 4 ko00944 Flavone and flavonol biosynthesis 0.0000184 4 ko00908 Zeatin biosynthesis 0.00016421 3 ko00062 Fatty acid elongation 0.00024383 7 ko00073 Cutin, suberine, and wax biosynthesis 0.00038521 6 ko00270 Cysteine and methionine metabolism 0.00047493 12 ko04075 Plant hormone signal transduction 0.00060961 16 ko00830 Retinol metabolism 0.00074745 3 ko00982 Drug metabolism-cytochrome P450 0.00093032 7 ko04540 Gap junction 0.00161368 6 ko00350 Tyrosine metabolism 0.00278399 5 ko01040 Biosynthesis of unsaturated fatty acids 0.00278399 5 ko00592 Alpha-linolenic acid metabolism 0.00609769 5 ko00980 Metabolism of xenobiotics by cytochrome P450 0.00609769 5 ko00140 Steroid hormone biosynthesis 0.00791195 3 ko00564 Glycerophospholipid metabolism 0.01020997 7 7 BioMed Research International    Table 2). The high expression quality (ExN50) of assembled contigs (N50) was revealed by majority of the contigs (>1,150 bp) ( Table 2).

Gene Expression and Differential Expression Analysis.
The overall gene expression was higher in L-species than in viviparous X-species observed by FPKM values (Figure 3(a)). The principal component analysis (PCA) indicated the close relation of samples within species while a relative high distance between samples from the two species (Figure 3(b)). It was supported by the results of the average Pearson's coefficient of correlation (Figure 3(c)), which indicates extensive genetic dissimilarities between the two species at the different developmental stages.
The expressed genes were further screened for their differential expression (DEGs) using DESeq2 analysis based on jlog 2 foldchangej ≥ 1 and false discovery rate ðFDRÞ < 0:05. Relatively higher numbers of DEGs were regulated in later growth stages than the early stages (Figures 4(a) and 4(b)). The differential expression analysis between the two species showed 5,559, 5,164, 6,375, and 7,391 DEGs for L1 vs. X1, L2 vs. X2, L3 vs. X3, and L4 vs. X4, respectively, resulting in a total of 10,956 unique DEGs among the total expressed unigenes for all four stages (Figure 4(b); Additional Table S1). The majority of genes were downregulated with the relative growth in next developmental stages from X1 to X4 (Figure 4(c)). By comparing all of the DEGs from all developmental stages, we identified 2,551 (23.3%) core conserved genes constantly differentially expressed between the two species (Figure 4(d); Additional Table S1). Along with core conserved DEGs, a high number of specific DEGs (20.7%, 2,267) were observed between X-and Lspecies at the 4 th developmental stage, correlating with the

11
BioMed Research International appearance of viviparous leaves in X-species. These genes may represent key genes involved in viviparity trait in water lily.

Functional Annotation and Enrichment Analysis of
DEGs. The differentially expressed genes were mapped to Gene Ontology (GO) terms in the GO database [25] for better understanding of functions and annotations in different developmental stages of viviparous and control Nymphaea species. GO functional enrichment analysis was performed adjusting P value of 0.05 as the cutoff (Additional Figure 1). A total of 4,943 GO terms were annotated to the 19,018 unigene hits (Additional Table S2). Among these terms, the maximum 52.26% (2,583) GO terms were belonged to the class "biological processes" (BP) followed by molecular functions (MF) (32.86%, 1,624 terms) and "cellular components" (CC) (11.77%, 582 terms). In CC, the most enriched GO terms were "nucleus" with 4,736 genes and membranous components including "integral components of membrane" with 3,909 genes and "plasma membrane" with 2,457 genes (Table 3), while in MF, most enriched terms were "ATP binding" with 3,123 genes followed by "metal ion binding" and "DNA binding transcription factor activity" with 2,793 and 1,062 genes, respectively. Among the biological processes, the "DNA integration" with 766 genes and "DNA recombination" with 735 genes were on top hit.
The DEGs were further evaluated for their functional enrichment between pairwise comparisons based on KEGG database [24]. The total 1,967 DEGs in the four developmental stages could be enriched in 173 unique KEGG pathways. At the early stages of plant development (L1 vs. X1 and L2 vs. X2), the DEGs related to "plant hormone signal transduction" and "flavonoid biosynthesis" pathways were downregulated, while the genes involved in "photosynthesis-antenna protein" and carbon fixation-  related pathways were significantly upregulated among L-and X-species (Additional figure 2). The DEGs related to "base excision repair," "DNA replication," and various metabolism pathways were downregulated at the third developmental stage of both L-and X-species species, while the DEGs for photosynthesis, carbon fixation, and other biosynthesis pathways were still upregulated (Additional figure 2). During the fourth developmental stage, the DEGs for "fatty acid elongation," "DNA replication," "cell cycle," and "meiosis" were differentially expressed (upregulated), while the genes for "flavonoid biosynthesis," "fatty acid metabolism," and various biosynthesis and metabolism-related pathways were also differentially expressed (downregulated) between the two species (Additional figure 1). Among all 173 annotated pathways, 127 pathways were conserved during all four developmental stages of X-and L-species (Table 4 and  Additional Table S3).

Plant Hormone Signal
Transduction. Among the highly conserved pathways in all four developmental stages, the plant hormone signal transduction pathway showed the highest differential regulation. Similar results have been reported in response to biotic stress in other plants [30]. The auxininduced proteins (AUX/IAA) were differentially regulated in X-and L-species at various stages. The auxin response factor (ARF) transcriptional factor was significantly upregulated

13
BioMed Research International which resulted in suppression of small auxin-up RNA (SAUR) gene family in L-species. The glycoside hydrolase 3 (GH3) expression was continuously decreased in L-species at later growth stages, indicating that L-species was unable to continue cell enlargement and plant growth ( Figure 5). Along with auxins, the cytokinin signaling-related gene histidine kinase (CRE1) was downregulated, and type A ARR was differentially expressed for signaling regulation. The downregulation of these genes in L-species is indicative of lower cytokinin signaling in L-species for cell division and shoot initiation in later developmental stages ( Figure 6). We further identified the differential regulation of cyclin-D3 (CYCD3) protein which is responsible for cell division and mitotic cycles in leaves [31]. The gene for MYC2 family was also downregulated in L-species which are well known to shape the plant growth and development [32]. Together with the higher auxin, cytokinin, and jasmonic acid contents in the leaves of X-species, the upregulation of key DEGs benefits the development of leaf outgrowth at later plant developmental stages. The genes related to Natriuretic Peptide Receptor 1 (NPR1) and pathogenesis-related protein (PR1) were upregulated in Lspecies which may be involved in the adenosine triphosphate (ATP) synthesis and responses to various stresses [33,34] ( Figure 5).

Carbon Fixation and Photosynthesis.
The plant development is also affected by the energy-related pathways including the changes in photosynthetic efficiency in response to modifications in photosynthesis-antenna proteins, carbon fixation pathways, and the photosynthesis pathway [35].
The modifications in photosynthesis and carbon fixationrelated pathways are accompanied by other energy and metabolism-related pathways including galactose metabolism, starch and sucrose metabolism, nitrogen metabolism, pentose phosphate pathway, citrate cycle, and carbon fixation in photosynthetic organisms [36]. The DEGs between L-and X-species were enriched in these pathways signifying large-scale transcriptional changes in energy-related pathways. Three DEGs in light harvesting chlorophyll protein complex (LHC) were observed to be significantly upregulated. Seven photosystem II (PSII) proteins (PsbK, PsbH, PsbI, PsbW, PsbZ, PsbO, and PbsQ) and one photosystem I (PSI) protein (PsaE) were differentially regulated. Only two genes related to "delta" and "a" F-type ATPase were downregulated. Four genes in photosynthesis electron chain transport (PetJ and PetF) were differentially expressed in the two species ( Figure 6).
Twenty DEGs (12 upregulated, 7 downregulated, and one up-/downregulated in at least one growth stage of development in the two species) were enriched in carbon fixation pathway in photosynthetic organisms ( Figure 6). Two dehydrogenase, four malate dehydrogenase (MDH), one NADP+ malate dehydrogenase genes were highly downregulated in L-species; however, other MDHs were upregulated. One FBPase and one SBPase genes were downregulated. Two aldose-related genes were downregulated at the 4 th developmental stage, while one phosphoenolpyruvate carboxylase (PEPC) gene was upregulated and one downregulated at the 4 th stage of plant development. Meanwhile, two phosphoribulose epimerase carboxylase (RuBPC) genes were    14 BioMed Research International downregulated, and one gene was upregulated throughout the plant development. The differential regulation patterns of various genes with the same annotation within the energy-related pathways showed complex differential transcription signatures in Nymphaea species. These expression changes in these pathways indicate their importance in the viviparous growth at later plant growth stages (Figure 7).

Phenylpropanoid/Flavonoid
Biosynthesis. Phenylpropanoids and flavonoids play vital roles in plant development by acting as essential components of cell walls, protectants against high light and UV radiation, phytoalexins against herbivores and pathogens, and floral pigments to mediate plant-pollinator interactions [37]. Nineteen DEGs were upregulated, eleven downregulated, and one up-/downregulated in the two Nymphaea species. Among the enriched proteins in these pathways, four beta-glucosidase, three chalcone synthase (CHS), one flavonol synthase (FLS), and one shikimate-hydroxycinnamoyl-transferase genes were downregulated. Two phenylalanine ammonia lyase (PAL), eight peroxidase, and four caffeate-O-methyltransferase genes were upregulated ( Figure 7). These complex mechanisms of transcriptome expression reveal the importance of these DEGs in the differential development of viviparity in X-species than in L-species.

Validation of RNA Analysis by qRT-PCR.
To further validate the expression of the identified DEGs between the two Nymphaea species at different developmental stages, we selected top ten DEGs (top five positively and top five negatively expressed genes) and performed qRT-PCR expression profiling. The qRT-PCR results of the selected genes were almost consistent with that of RNA-seq analysis ( Figure 8). There was significant difference (P ≤ 0:05) of expression levels between the two studied species for all evaluated genes at the four growth stages. This result supports the DEG analysis and subsequent interpretations.

Transcriptome Enrichment and Genetic Basis of
Viviparity in Water Lily. The extrachromosomal genome [38], transcriptome, and proteome [10] analyses for various morphological and physiological traits have been reported in water lilies. However, this is the first research focused on the molecular basis of viviparous growth in water lilies. We generated the whole transcriptome sequence from leaf tissues at various plant developmental stages, assembled the sequencing data and annotated the differentially expressed for functional assessment. The vivipary reproduction has a great flexibility for parental control of embryonic development, which in turn allows viviparous organisms to reproduce successfully in adverse environments [39]. There are some studies reported for the genetic bases of transition from oviparity to viviparity in animals by comparing genomic and transcriptomic data [39], but our knowledge of viviparity trait in plants is still limited. We used the N. micrantha as a model viviparous water lily to compare its transcriptome with the nonviviparous N. colorata at various developmental stages. The highest number of DEGs between the two species was observed at the 4 th developmental stage of water lily. The viviparity phenomenon was observed to be associated with 127 highly conserved metabolic pathways in the four developmental stages. The top enriched pathways belong to plant hormone signal transduction, DNA replication, cell cycle, photosynthesis, and carbon fixation. The expression analyses of differentially expressed genes (DEG) indicate the complex network underlying the viviparous growth [40]. Besides, there were 3,574 DEGs conserved between the third and fourth growth stages that may be involved in plantlet development and may provide clues for further studies on identifying inducible/specific mechanisms/components involved in adaptation of viviparous growth.

Effect of Phytohormones on Viviparous Growth in Water
Lily. As plants grow, they develop new organs as primary and secondary leaves, lateral roots, and flowers [41]. The plant development is regulated by complex hormone interaction and signaling which helps various species to evolve in the wide range of environmental conditions [41]. Various biotic/abiotic stress combinations develop a new type of signal and response in plants, resulting in a novel transcription signature [42]. Through our comparative transcriptome study among the viviparous and nonviviparous water lily species, we found similar results as the biotic stress responses [43]. The biotic/abiotic stress can cause the variation in phytohormonal balance in plants. We identified the significant differential regulation of AUX/IAA, ARF, GH3, and SAUR. Auxin/Indole-3-Acetic acid (Aux/IAA) genes are the early response genes that trigger gene reprogramming precisely and rapidly under stress [44]. Auxin response factor (ARF) transcription factors are activated upon auxin perception and initiate downstream signaling pathways including the small auxin upregulated RNA (SAUR) genes [41]. SAURs regulate many auxin-mediated responses, specifically the tissue growth via cell elongation [41]. We identified the highly conserved genes for ARF (NC0286500) and SAURs (NC14G0174900, NC263970, NC7G0235790, and NC9G0168130) ( Figure 5). These genes may directly be involved in the development of viviparous outgrowth in water lilies by their involvement in cell division, enlargement, and differentiation [44]. The auxin in association with cytokinins also helps in cell differentiation [41]. The higher expression of CRE cytokinin may result in viviparous cell differentiations as CRE is known for acting on shoot apical meristem differentiation [45]. The changes in other signaling hormones also have their indirect role to enhance the viviparous plantlet growth. It is possible because CYCD3 cyclins was previously reported for mitotic cell divisions [31]. The upregulated expression of DEGs related to jasmonic acid (JA) signaling, i.e., MYC2, is also an indicator of high response to light phytochromes [46] that was reduced in Lspecies. It can finally be concluded that the viviparous plantlet differentiation could be due to the higher regulation of auxin and cytokinin and subsequent activation of other relevant phytohormones ( Figure 5). In previous studies, the role of abscisic acid (ABA) and gibberellic acid (GA) has also been reported for development of viviparity and influence 15 BioMed Research International on meristematic tissue, but in the present study, no DEG was observed to be involved in ABA or GA-based function for viviparity.

Involvement of Viviparous Plantlet in Photosynthesis.
Development of viviparous plantlet growth significantly affects the photosynthesis and subsequent energy-related processes [47]. The increasing expression of light harvesting chlorophyll protein complex (LHC3) and LHC5 could be a subsequent regulation effect of viviparity. Further, the downregulation of relatively higher number of photosynthesis associated proteins (photosystem II (PSII) and photosystem I (PSI) proteins) in L-species reduced its photosynthetic potential that possibly affected downstream energy metabolism-related pathways [48,49]. The PsbQ and PsbO proteins of PSII-complex were downregulated in L-species, while they have higher expression in X-species ( Figure 6). These proteins have been reported for their role to stabilize the interaction between the membrane-bounded PSII subunits and the related proteins, i.e., PsbP. Hence, they may be involved in the regulation and evolution of photosynthesis in viviparous plantlet.

Role of Flavonoid Pathways in Viviparous Water Lily.
Phenylpropane and peroxidases are the precursor of lignin biosynthesis, which are an essential element of cell wall [50]. The phenylpropanoids are a group of plant secondary metabolites derived from phenylalanine and have a wide variety of functions both as structural and signaling molecules [51]. Lignin is derived from phenylalanine and has a wide variety of functions both as structural and signaling molecules [51,52]. The DEGs for peroxidases, caffeate Omethyltransferase, CYP75B1, and PAL were upregulated in the nonviviparous species (Figure 7) of water lily. These secondary metabolites are known to play key roles in inhibition of seed germination and reduced water permeability [50]. Inversely, the structural genes (CHS) [53] and genes related to beta glucosidase and FLS which are known for tissue development and the anthocyanin biosynthesis [54,55] showed downregulation in nonviviparous species (Figure 7). The higher expression of these genes in viviparous species revealed their importance in the viviparity phenomenon.

Conclusions
The comparison of transcriptome of viviparous species N. micrantha and nonviviparous species N. colorata of water lilies revealed the variation in expression of various genes indicating their putative role in viviparity. Nymphaea micrantha showed upregulated genes for plant hormone signal transduction including AUX/IAA, ARF, GH3, and SAUR gene families. Regulation of these genes involved in cell division, elongation, and differentiation showed their association to viviparity in N. micrantha. The increased expression levels of these genes in N. micrantha triggered downstream phytohormone signaling cascade as noticed by the regulation of genes such as MYC2s and CYCD3. The viviparity in plants also modulates changes in the important energy-related pathways. Genes related to both light harvesting chlorophyll complex and photosystem I and photosystem II were differentially expressed in both species. Importantly, increased expression of major genes related to the above cited pathways in N. micrantha possibly regulates the phenomenon of viviparity. Overall, the transcriptomic signatures identified in this study are important basis for future research of viviparity in water lilies and other plant species.

Data Availability
The RNA-seq raw data has been submitted to NCBI GEO under the accession number GSE164888.

Disclosure
The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Supplementary Materials
Additional Figure 1: list of top enriched Gene Ontology (GO) terms between two Nymphaea species N. micrantha and N. colorata, at four developmental stages. Additional Figure 2: list of top twenty enriched KEGG pathways between two Nymphaea species N. micrantha and N. colorata, at four developmental stages. Additional Table S1: list of unique differentially expressed genes among four developmental stages of N. micrantha (X) and N. colorata (L) species. Available at https://figshare.com/articles/dataset/ Additional_information_tables_xlsx/19688034). Additional Table S2: list of unique Gene Ontology (GO) terms annotated to the differentially expressed genes among four developmental stages of N. micrantha (X) and N. colorata (L) species. Available at https://figshare.com/articles/dataset/ Additional_information_tables_xlsx/19688034. Additional 16 BioMed Research International